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PREFACE 


This  technical  report  was  prepared  by  the  Research  Staff  of  the  Electrical 
Engineering  Department,  School  of  Engineering,  the  The  University  of  Alabama 
in  Huntsville.  The  purpose  of  the  report  is  to  provide  documentation  of  the 
technical  work  performed  and  results  obtained  under  delivery  order  0007  of 
MICOM  Contract  No.  DAAH01-82-D-A008 ;  Dr.  N.  A.  Kheir,  Principal  Investigator. 

The  project  documented  herein  was  performed  by  Patrick  A.  Tilley  and 
Donn  Hall.  Dr.  M.  M.  Hallum,  III,  Chief,  Systems  Evaluation  Branch,  Army 
Missile  Laboratory,  US  Army  Missile  Command,  was  technical  monitor,  and 
Mr.  D.  L.  Cobb  provided  technical  coordination. 

The  authors  wish  to  acknowledge  the  valuable  discussions  and  assistance 
provided  throughout  the  task  by  D.  L.  Cobb  and  Dr.  M.  M.  Hallum  of  the 
Systems  Evaluation  Branch. 

The  2.75  Rocket  6-DOF  simulation  program,  described  in  Section  IV,  was 
developed  by  H.  Lanier,  MICOM,  and  modified  by  R.  A.  Dillard,  MICOM. 

The  technical  viewpoints,  opinions  and  conclusions  expressed  in  this 
report  are  those  of  the  authors  and  do  not  necessarily  express  or  imply 
policies  or  positions  of  the  US  Army  Missile  Command. 


I.  INTRODUCTION 

The  primary  goal  of  this  task  was  to  find  impact  points  of  the  second 
stage  and  adapter  of  the  PERSHING  II  (PXI)  missile.  This  task  was  necessary 
because  an  accu^te  approximation  of  major  debris  impact  points  was  needed 
for  PERSHING  II  test  flights.  To  accomplish  this  goal,  the  MSIXB1  6-DOF 
Simulation  was  the  primary  tool  used.  It  is  described  in  Section  II.  The 
results  of  debris  impact  studies,  using  MSIXB1,  are  then  presented  in  Section 
III.  Additionally,  another  method  of  debris  evaluation,  the  2.75  Rocket  6-DOF 
simulation,  is  described  in  Section  IV  to  facilitate  future  use  in  debris 
studies.  Then,  in  Section  V,  applicable  plotting  programs  are  described. 

II.  DEBRIS  SIMULATION  PROGRAM  (MSIXB1) 

A.  Introduction 

The  method  chosen  to  approximate  impact  points  uses  MSIXB1  and  input 
missile  characteristics  form  the  U-70  6-DOF  Simulation  of  the  PII  missile. 
Specifically,  an  input  file  containing  values  such  as  velocity  and  position 
is  used  by  MSIXB1  to  find  an  impact  point  for  each  set  of  inputs.  Thrust 
tables  generated  for  various  cutoff  times  are  also  used  by  the  program. 

The  main  MSIXBl  program  tracks  the  second  stage  to  impact  while  an 
alternate  path  in  the  program  can  track  the  adapter  or  any  other  debris. 

MSIXBl  finds  the  impact  points  by  thrusting  until  cutoff  and  then  letting  the 
debris  fall  ballistically  until  impact.  Inputs  can  be  varied  to  supply  a 
footprint  of  impact  points  which  can  be  analyzed  statistically. 

B.  Data  Inputs 

Data  changes  are  input  into  MSIXBl  in  an  efficient  manner.  The  type 
of  data  is  determined  by  the  first  character  In  each  line  of  data.  If  this 
character  is  numeric,  the  line  of  data  is  interpreted  as  input  data  to  be 
stored  in  an  array  which  will  be  used  by  the  program.  If  the  first  character 
is  alphabetic,  it  will  trigger  the  program  to  run,  to  go  into  its  ballistic 
phase,  or  to  end  execution. 

After  a  line  of  data  is  detected  an  input  data,  the  first  number  read 
defines  the  element  number  of  the  input  array.  Input  data  is  begun  a^ter  a 
space  is  encountered.  Two  methods  can  then  be  used  to  read  the  data  dn  the 
line.  The  first  method  records  the  second  number  as  the  value  of  the  element 
defined  by  the  first  number.  All  of  the  following  numbers  are  then  entered 
into  the  data  array  sequentially.  An  example  is  given  below  with  an  explana¬ 
tion  presented  in  Table  I. 

73  15.2  16.3  17.4  3.277  1086.0  -107.32  0. 
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TABLE  I.  MSIXB1  SINGLE  VARIABLE  INPUT  SCHEME  EXAMPLE 


V 

•  i 

k 


1,1 


r. 

>, 


Data  Array  Element  Number 

Value 

73 

15.2 

74 

16.3 

75 

17.4 

76 

3.277 

77 

1086.0 

78 

-107.32 

79 

0. 

The  second  method  allows  a  table  of  values  to  be  read  in  from  one  or 
more  lines  in  the  input  data  file.  The  first  number  is  again  used  as  the 
beginning  data  array  element  number.  The  second  number  indicates  the  number 
of  elements  in  the  table  of  values,  while  the  third  number  is  set  to  zero  in 
order  to  align  the  inputs  correctly.  The  fourth  number  is  the  first  number 
in  the  table  of  values  and  is  placed  in  the  data  array  as  an  element  defined 
by  the  first  number.  The  remaining  values  are  then  entered  into  the  data 
array  in  sequential  order  as  above.  An  example  follows  with  the  explanation 
given  in  Table  II. 

400  14  0  5.  77.32  10.  89.15  15.  97.98  20.  123.50  25.  141.34  30.  136.2  35.  -17.1 
TABLE  II.  MSIXB1  TABLE  INPUT  SCHEME  EXAMPLE 


Data  Array  Element  Number 

Corresponding  Array  Value 

400 

5. 

401 

77.32 

402 

10. 

403 

89.15 

404 

15. 

405 

97.98 

406 

20. 

407 

123.50 

408 

25. 

409 

141.34 

410 

30. 

411 

136.2 

412 

35. 

413 

-17.1 

a 
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The  first  number  in  each  set  of  active  data  represents  a  time  corres¬ 
ponding  to  the  second  number.  Tables  are  called  in  an  interpolation  routine, 
and  an  interpolated  value  is  calculated  according  to  the  surrounding  times 
and  their  values. 

The  alphabetic  data  inputs  mentioned  earlier  are  available  to  perform 
special  operations  in  the  program.  Presently,  only  three  alphabetic  inputs 
are  tapped  for  use  in  the  program.  These  are  "R"  for  running,  "C"  for  going 
into  ballistics,  and  "E"  for  ending.  Operation  of  the  program  ceases  only 
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when  an  "E"  is  encountered.  Thus,  the  data  array  can  be  changed  at  any  time, 
and  the  program  can  be  rerun  using  only  one  input  data  file.  An  example  is 
given  in  Table  III. 

TABLE  III.  MSIXB1  ALPHABETIC  INPUT  SCHEME  EXAMPLE 


LO.  12. 


161  7.  15.  32. 
161  10.  25. 

R 


data  array  values 

program  runs  using  above  values 

program  goes  into  ballistics  after  run  is 

completed 

new  data  array  values 

data  array  values  replacing  part  of  above  line 
program  runs  using  new  values  161=10. 

162=25.  163=32. 

program  ends. 


As  shown  in  the  above  example,  only  the  final  values  corresponding  to  data 
array  elements  which  precede  each  "R"  are  used  during  the  run.  This  proves 
to  be  a  very  useful  tool  since  an  existing  input  data  file  can  be  altered  by 
simply  adding  a  new  line  of  data  after  the  old  line.  This  is  convenient  since 
experimentation  with  new  inputs  can  be  made  without  disturbing  functioning 
inputs.  A  list  of  all  possible  active  inputs  is  given  in  Tables  IV-VI.  The 
English  system  is  used,  and  time  inputs  have  units  of  seconds. 

C.  Aerodynamics 

MSIXB1  utilizes  aerodynamics  based  on  coordinate  transformations 
using  quaternions.  The  four  quaternion  parameters  are  calculated  in  the  sub¬ 
routine  KINE.  Quaternions  are  an  analogy  to  complex  numbers,  but  use  three 
components  instead  of  two.  They  are  used  because  of  their  excellent  co¬ 
ordinate  system  rotation  properties  which  reduce  the  number  of  calculations 
required.  QAP(l),  QAP(2),  QAP(3),  and  QAP(4)  are  quaternions  used  to  find 
coefficients  for  coordinate  transformations.  These  coordinate  transformation 
coefficients  are  then  used  in  subroutine  AERO  to  properly  simulate  the  missile 
motion. 

The  coordinate  system  used  in  MSIXB1  is  shown  in  Figure  1.  Also  shown 
in  Figure  1  is  the  U70  coordinate  system,  the  other  coordinate  system  used  in 
the  PERSHING  II  study.  The  x-axis  is  downrange  in  both  systems,  and  the 
following  transformations  are  needed. 


> _■*  V  V  *.1 


TABLE  IV.  MSIXB1  SINGLE  INPUTS 


Array  Element  Number 

Variable 

Subroutine 

1 

X  Position  (Ft) 

KINE 

2 

Y  Position  (Ft) 

KINE 

3 

-Z  Position  (Ft) 

KINE 

4 

Total  Velocity  (Ft/S) 

KINE 

5 

Initial  Yaw  Angle  (Deg) 

KINE 

6 

Initial  Pitch  Angle  (Deg) 

KINE 

7 

PSI  -  Yaw  Euler  Angle  (Deg) 

KINE 

8 

THETA  -  Pitch  Euler  Angle  (Deg) 

KINE 

9 

PHI  -  Roll  Euler  Angle  (Deg) 

KINE 

10 

Roll  Rate  (Deg/S) 

KINE 

11 

Yaw  Rate  (Deg/S) 

KINE 

12 

Pitch  Rate  (Deg/S) 

KINE 

13 

pi stance  Between  Force  &  CG  in  Z 
Direction  (In.) 

AERO 

14 

Distance  Between  Force  &  CG  in  Y 
Direction  (In.) 

AERO 

15 

Acceleration  of  Gravity  (Ft/S2) 

MAIN,  AERO 

19 

Weight  (Lbs) 

AERO 

20 

CG  Mass  Property  Table  Number 

AERO 

21 

X  Mass  Property  Table  Number 

AERO 

22 

Y  Mass  Property  Table  Number 

AERO 

23 

Moment  arm  used  to  calculate 
Variable  M2 

AERO 

25 

Print  Option  for  Metrics,  Etc. 

MAIN 

26 

Total  Body  Rate  Limit  0  =  No  Limit 

MAIN 

27 

Print  Option  for  Delta  Values 

MAIN 

28 

Roll  Rate  Limit 

MAIN 

29 

QR  Rate  Limit 

MAIN 

30 

Iteration  Time 

MAIN 

31 

Print  Iteration  Interval 

MAIN 

32 

Alternate  Print  Iteration  for 
TIME.GT.  A( 205 ) 

MAIN 

33 

Time  Limit  For  Run 

MAIN 

34 

Seconds  After  Separation  That 
Thrust  Reversers  Are  Activated 

AERO 

38 

Distance  Between  Center  of  Gravity 
And  Thrust  Force  to  Calculate 
Variable  M3 

AERO 

40 

Print  Option  for  Deltas 

MAIN 

41 

Variable  Delta  Option  Time 

MAIN 

42 

X  Velocity  Vector  Added  to  Debris 
After  Time  A (49) 

MAIN 

43 

Y  Velocity  Vector  Added  to  Debris 
After  Time  A (49) 

MAIN 

44 

Z  Velocity  Yector  Added  to  Debris 
After  Time  A(49) 

MAIN 
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TABLE  IV  (cont'd).  MSIXB1  SINGLE  INPUTS 


— 

Array  Element  Number 

Variable 

Subroutine 

45 

Ballistic  Coefficient  of  Debris 

BALLS! 

46 

Iteration  Time  For  Ballistic  Subroutine  BALLS!] 

47 

Print  Iteration  For  Ballistic  Subroutine  BALLST 

49 

Time  After  Separtion  When  Secondary 

Debris  is  Tracked 

MAIN 

50 

IN-LB/PSI  Coefficient  For  Roll 

AERO 

51 

IN-LB/PSI  Coefficient  For  Pitch-Yaw 

AERO 

52 

Pitch  Yaw  Phi  Angle 

AERO 

90 

Print  Start  Time 

MAIN 

101 

X  Offset  Distance  Between  A ( 1 ) 

&  Start  Point 

MAIN 

102 

Y  Offset  Distance  Between  A (2 ) 

&  Liftoff  Point 

MAIN 

160 

Weight  From  Which  Thrust  Ratio 

is  Calculated 

AERO 

161 

Moment  1  Input 

AERO 

162 

Moment  2  Input 

AERO 

163 

Moment  3  Input 

AERO 

164 

Upper  Time  Limit  At  Which  Moment 

2  Equals  0 

AERO 

165 

Lower  Time  Limit  At  Which  Moment 

2  Equal s  0 

AERO 

205 

Time  Which  Changes  Print  Iteration 

to  A ( 32 ) 

MAIN 

•  777^  '."v  777  ^  <•. 
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TABLE  V.  MSIXB1  INPUTS  IN  AERO  >: 


, 


Table  Array  Start  Number 

Table  Subject 

*i 

305 

Delta  Pitch  Attack  Angles 

\*V 

320 

Delta  Yaw  Attack  Angles 

500 

Forward  Thrust  Before  TR  (Thrust  Reversal) 

'•A 

550 

Reverse  Thrust  Before  TR 

600 

Rate  of  Weight  Change  Before  TR 

•J 

vi 

1500 

Forward  Thrust  After  TR 

«.■) 

1550 

Reverse  Thrust  After  TR 

*• 

1600 

Rate  of  Weight  Change  After 

TR 

1700 

Thrust  of  Dome  Reverser  Ports 

nJ 

2600 

CG  Mass  Property 

1-4 

■1 

* 

2620 

X  Mass  Property 

2640 

Y  Mass  Property 

TABLE  VI.  MSIXB1  ALPHABETIC  INPUTS 

£ 

% 

Letter 

Command 

C 

Calls  BALLST  After  Rur 

1 

E 

Ends  Run 

r'l 

R 

Runs  Program 

r. 

’  • 

« 

*  « 

•J 

•J 
• , 

u 

n 

%  1 

0 

•: 

'1 

y 


g 
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« 
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.* 
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x  =  x,  y  -  z,  y  =  -y 


z 


Figure  1.  U70  and  MSIXB1  coordinate  systems. 

However,  it  should  be  noted  that  the  negative  sign  for  z  is  taken  into  account 
when  the  z  value  is  entered  into  the  program  as  -A(3) . 

D.  Ballistic  Phase 

New  subroutines,  BALLST  and  ATMD,  were  developed  and  added  to  the  basic 
MSIX  program  to  permit  simulation  of  missile  debris.  The  only  forces  acting 
upon  debris  are  gravity  and  drag.  The  weight  of  the  debris  is  taken  into 
account  by  the  ballistic  coefficient  "BC,"  which  is  weight  divided  by  drag 
coefficient.  Since  the  debris  begins  falling  from  a  point  outside  the  atmos¬ 
phere,  ATMD  has  been  included  to  supply  the  changing  density  of  the  air  which 
affects  drag.  Time  and  print  iteration  intervals,  specifically  for  the  ballis¬ 
tic  phase,  are  included  as  inputs  A(46)  and  A(47) .  Iterations  continue  until 
the  z  component  of  the  debris'  position  is  greater  than  zero.  This  impact 
point  is  then  printed  in  both  English  and  metric  coordinates,  and  the  run  is 
terminated. 

E.  Program  Operation 

A  complete  set  of  MSIXB1  files  consists  of  a  CSS  file,  a  FTN  (fortran) 
file,  and  TSK  (task)  file.  Alterations  to  the  program  are  made  in  the  FTN 
file.  This  FTN  file  is  compiled  into  the  machine  code  which  comprises  the 
TSK  file.  The  CSS  file  is  a  file  which  contains  instructions  to  run  the  TSK 
file.  This  CSS  file  is  initiated  by  simply  typing  MSIXB1,  the  name  of  data 
file  to  be  used,  and  return.  The  CSS  program  is  presented  in  Figure  2. 


L  MSIXB1 
AS  1,C0N: 

AS  3, NULL: 

AS  5.01.DAT 
AS  6, CON: 

AS  7, CON: 

AS  8, CON: 

AS  9, CON: 

AS  4, CON: 

ST 

$EXIT 

Figure  2.  MSIXB1  CSS  file 
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The  data  file  name  should  have  the  extension  .DAT  to  function  properly.  A 
flow  chart  of  the  FTN  file  is  given  in  Appendix  A,  and  the  FTN  file  itself 
is  presented  in  Appendix  B. 

III.  DEBRIS  IMPACT  STUDIES 

A.  Introduction 

Three  pieces  of  debris  have  been  studied  thus  far.  These  are  the 
second  stage,  the  adapter  between  the  second  stage  and  the  re-entry  vehicle, 
and  the  dome  closures  covering  the  thrust  reversal  stacks.  The  studies  used 
several  different  separation  times  for  each  piece  of  debris.  Thus,  several 
sets  of  initial  conditions  had  to  be  used  in  the  data  files.  Also,  various 
initial  velocity  vectors  were  added  to  yield  a  circle  of  impact  data.  Input 
A(49)  controls  the  time  at  which  the  secondary  debris  begins  falling.  The 
value  that  is  used  to  describe  the  weight  and  drag  properties  of  the  object 
is  the  ballistic  coefficient.  This  value  is  defines  as  weight  divided  by  the 
drag  coefficient. 

B.  Adapter  Study 

A  ballistic  coefficient  of  8.8  was  used  for  the  adapter.  Four  sets 
of  added  initial  velocity  vectors  were  input,  and  a  run  was  made  for  each  set 
at  three  different  separation  times.  The  vectors  were  all  combinations  of 
positive  and  negative  vectors  added  in  z  and  y  directions.  Several  different 
impulse  moments  were  used  on  the  missile  so  that  a  circular  error  pattern  could 
be  found  for  each  set  along  with  an  impact  footprint.  The  adapter  impacted  at 

shorter  range  and  had  a  smaller  cross-range  deviation  than  did  the  second 

stage.  Thus,  the  debris  stopped  short  in  downrange  causing  a  smaller  deviation 
in  the  y-direction. 

C.  Port  Closure  Study 

Two  different  weight  port  closures  were  studied  which  had  ballistic 
coefficients  of  7.4  and  22.2.  Chamber  prelsure  of  the  second  stage  was 
analyzed  to  yield  the  average  and  limiting  values  of  added  initial  velocity 

vectors.  Twelve  runs  were  made  for  each  of  eight  separation  times  to  fully 

evaluate  the  debris  pattern  of  each  piece  of  debris.  Plots  were  made  showing 
the  characteristic  impact  ellipse  for  each  debris  mass,  with  the  boundary  of 
the  ellipse  being  determined  by  the  limiting  initial  added  velocities. 

D.  Second  Stage  Impact  Study 

The  second  stage  study  was  similar  to  the  secondary  debris  studies 
except  that  forces  continued  to  act  on  the  body  even  after  thrust  cutoff. 

These  thrusts  occur  due  to  the  release  of  pressure  through  the  thrust  reversal 
ports.  The  program  simulated  motion  until  all  thrusts  were  zero.  Then 
ballistic  and  atmospheric  models  permit  simulation  from  this  cutoff  point 
until  impact. 

Runs  using  fifty  random  moments  were  made  for  the  second  stage  with  various 
velocity  vectors  added.  A  circular  error  pattern  routine,  which  found  99% 
and  50%  error  circles,  was  used  to  analyze  this  data. 


IV.  2.75  ROCKET  6-DOF  SIMULATION  PROGRAM 


A.  Introduction 

The  2.75  Rocket  is  a  multi-mode  rocket.  It  can  function  as  either 
a  surface-to-surface  or  air-to-surface  rocket.  Additionally,  it  can  be  fired 
from  either  fixed  wing  aircraft  or  from  helicopters.  Its  simulation  is  based 
on  the  most  difficult  of  the  three  launch  platforms,  the  helicopter.  The 
simulation  is  versatile  enough,  though,  to  easily  adapt  to  the  other  plat¬ 
forms.  And  this  versatility  also  makes  the  simulation  useful  for  debris 
study. 


B.  Algorithm 

The  algorithm  for  the  simulation  is  common  although  some  of  the  pro¬ 
cesses  are  uncommon.  A  flowchart  and  complete  program  listing  can  be  found  in 
Appendix  C  and  D,  respectively. 

The  first  step  of  the  algorithm  is  to  predefine  the  values  of  all  varia¬ 
bles  at  launch.  A  "cycle"  begins  when  the  rocket  is  launched  and  starts  with 
an  initial  update  of  variables.  Added  to  these  variables  are  first  the  effect 
of  wind  and  rotor  downwash,  if  applicable,  and  then  the  corrections  for  aero¬ 
dynamic  and  atmospheric  conditions.  From  these  values  and  the  previous  values 
derivatives  of  the  variables  (Table  VII)  are  defined  and  calculated.  Finally, 
a  Runge-Kutta  integration  routine  is  used  to  integrate  the  derivatives 
(Table  VIII),  thereby  completing  the  cycle. 

C.  Coordinate  System 

There  are  three  coordinate  systems  to  be  considered  for  this  simula¬ 
tion.  These  are  centered  on  (1)  the  launch  platform,  (2)  the  rocket,  and 
(3)  the  ground.  All  three  have  the  same  conventions  and  are  shown  in  Figure  3 
It  should  be  noted  that  altitude  is  positive  downward. 


Figure  3.  2.75  uniform  coordinate  systems. 


Table  VII  is  a  list  of  the  derivative  array  and  the  derivatives' 
definitions.  Table  VIII  is  a  list  of  the  variables  that  are  the  result  of 
integrating  the  derivatives  array  and  their  definitions. 

TABLE  VII.  2.75  DERIVATIVE  ARRAY 


XX(1,1)  -  Angular  acceleration  about  the  missile  X-Axis 

XX( 1 ,2 )  -  Angular  acceleration  about  the  missile  Y-Axis 

XX(1,3)  -  Angular  acceleration  about  the  missile  Z-Axis 

XX (1,4)  -  Linear  acceleration  along  the  missile  X-Axis 

XX( 1 ,5 )  -  Linear  acceleration  along  the  missile  Y-Axis 

XX ( 1 ,6)  -  Linear  acceleration  along  the  missile  Z-Axis 

XX (1,7)  -  Linear  velocity  along  the  ground  X( range) -Axis 

XX{1,8)  -  Linear  velocity  along  the  ground  Y(deflection)-Axis 

XX(  1 ,9)  -  Linear  velocity  along  the  ground  Z(vertical)-  Axis 

XX ( 1,10)  -  Derivative  of  the  Euler  angle  PSI 

XX ( 1,11)  -  Derivative  of  the  Euler  angle  THETA 

XX( 1,12)  -  Derivative  of  the  Euler  angle  PHI 

XX(1,13)  -  Derivative  of  the  mass  (time  rate  of  change  of  mass) 

XX(1,14)  -  Acceleration  along  launcher  center  line 

XX(1,15)  -  Velocity  along  launcher  centerline _ 


TABLE  VIII.  2.75  STATE  VARIABLE  ARRAY 

XX(2,1)  -  Angular  velocity  about  missile  X-Axis  (Roll  Rate) 

XX(2,2)  -  Angular  velocity  about  missile  Y-Axis  (Pitch  Rate) 

XX (2, 3)  -  Angular  velocity  about  missile  Z-Axis  (Yaw  Rate) 

XX(2,4)  -  Linear  velocity  along  missile  X-Axis 

XX (2,5)  -  Linear  velocity  along  missile  Y-Axis 

XX(2,6)  -  Linear  velocity  along  missile  Z-Axis 

XX(2,7)  -  Position  with  respect  to  ground  X-Axis  (Range) 

XX(2,8)  -  Position  with  respect  to  ground  Y-Axis  (Deflection) 

XX(2,9)  -  Position  with  respect  to  ground  Z-Axis  (Altitude) 

XX(2,10)  -  Euler  angle  PSI 

XX (2,11)  -  Euler  angle  THETA 

XX (2,12)  -  Euler  angle  PHI 

XX( 2,13)  -  Missile  mass 

XX(2,I4)  -  Velocity  along  launcher  centerline 

XX( 2 ,15 )  -  Position  with  respect  to  launcher _ 

It  is  important  to  note  that  XX(1,1)  does  not  go  through  the  integration 
routine.  The  roll  rate  XX(2,1)  is  found  through  the  use  of  a  look-up  table 
based  on  the  time  since  launch. 


In  the  helicopter  system,  all  of  the  rocket's  launch  conditions  are 
calculated  from  the  position  of  the  helicopter.  Important  angles  in  the  heli¬ 
copter  system  are  illustrated  in  Figures  4  and  5.  More  information  on  the 
helicopter  as  a  launch  platform  can  be  obtained  from  Reference  [1]. 


Figure  4.  AH-1G  helicopter  angles  (side  view). 

The  numbered  angles  in  Figure  4  are  defined  as  follows: 

1.  Total  horizontal  angle  of  attack 

2.  Dive  angle 

3.  Helicopter  Altitude 

4.  Quadrant  elevation  of  launcher 

The  missile  coordinate  system,  illustrated  in  Figures  6  and  7,  is  the 
system  in  which  all  of  the  flight  characteristics  of  the  missile  and  all  of 
the  external  forces  on  the  missile  are  calculated.  Figure  5  illustrates  the 
rocket  with  respect  to  the  launcher  tube  shortly  after  launch.  Figure  7 
illustrates  the  basic  body  angles  and  velocity  vectors  of  the  rocket. 

The  ground  coordinate  system  (Figure  8)  is  based  on  the  position  of  the 
target.  Through  use  of  this  system,  the  inertial  position  of  the  rocket  is 
determined  as  well  as  the  position  of  impact.  The  simulation  assumes  a  flat 
earth  and  a  uniform  gravitational  field. 

Transformation  between  the  helicopter  and  missile  systems  is  performed 
simply  by  adding  and  subtracting  predefined  distances  based  on  parameters 
existing  at  the  time  of  calculation.  The  transformation  between  the  rocket 
and  the  earth-fixed  coordinate  systems  is  more  difficult,  and  the  process 
used  in  this  simulation  is  rather  unusual. 


DDIS 

XMW 

ZMW 

VXZ 

THEO 


Position  with  respect  to  laucher 
X  component  of  velocity  vector 
Z  component  of  velocity  vector 
Velocity  in  XZ  plane 
Total  missile  attitude  with 
respect  to  horizontal 


Figure  6.  2.75  Rocket  coordinate  system  with  respect  to  launcher. 


Figure  7.  2.75  Rocket  coordinate  system  body  angles  and  velocity  vectors. 


XMXH  =  X  component  of  inertial  position 
YMYH  =  Y  component  of  inertial  position 
ZMZH  =  Z  component  of  inertial  position 
VRF  =  Relative  velocity  vector 


EARTH 


Figure  8.  2.75  Rocket  earth  fixed  coordinate  system. 

Usually,  the  derivatives  of  the  Euler  angles  are  repeatedly  integrated  to 
update  the  coordinate  transformation  matrix  (Data  Conversion  Matrix-DCM)  in 
terms  of  sines  and  cosines  of  the  Euler  angles.  However,  this  simulation 
digressed  from  that  method.  Integration  of  the  Euler  angle  derivatives  is 
used  to  update  the  Euler  derivatives  via  an  altered  DCM  array.  Thus,  the 
body  rates  are  used  to  refresh  the  DCM  matrix  instead  of  the  basic  Euler 
angles,  eliminating  one  step  in  each  calculation. 

D.  Downwash  Profile  Description 

During  the  first  few  seconds  of  the  rocket's  flight,  after  launched 
from  a  helicopter,  several  transients  associated  with  the  launch  platform  be¬ 
come  involved.  These  transients  include  vibration,  translation  and  rotation 
of  launch  platform  and  rotor  downwash.  As  soon  as  the  rocket  clears  the 
launcher  tube,  the  first  two  of  these  transients  cease  to  be  involved.  The 
downwash  affects  the  rocket  until  the  time  that  the  rocket  clears  the  rotor 
blades,  and  the  effect  changes  with  respect  to  helicopter  velocity.  These 
effects  have  been  studied  at  MICOM  for  many  years  and  are  modeled  in  the  pro¬ 
gram  through  subroutine  HELVEL.  It  uses  a  data  pack  containing  the  downwash 
characteristics  of  an  AH-1G  helicopter.  Refer  to  Reference  [1]  for  more 
information  on  these  transients  and  their  effects. 

E.  Program  Variables 

The  purpose  of  this  section  is  to  define  the  three  types  of  variables 
in  the  2.75  digital  simulation.  They  are  input,  computed,  and  output  variables. 
Due  to  the  number  of  variables  and  the  uncertainty  of  some  of  their  definitions, 
there  is  a  possibility  that  a  complete  listing  is  not  presented  in  this  section. 
However,  all  efforts  have  been  made  to  ensure  that  the  definitions  given  are 
correct.  Any  variable  that  is  given  in  parenthesis  without  a  definition  is 
done  so  because  a  definition  is  not  known  at  this  time.  All  variable  lists 
are,  for  the  most  part,  in  the  order  that  they  appear  within  the  program. 


1.  Input  Variables 


Table  IX  is  a  list  of  the  variables  input  to  the  2.75  digital 
simulation  with  their  definitions.  (NAMELIST/PUTT/) 

TABLE  IX.  2.75  INPUT  VARIABLES 


- m us - 

ALTIN  -  Initial  missile  altitude 

DLETT  -  Delta  time  (step  size)  for  integration  routine 
DEL02  -  Delay  time 

DIVANG  -  Angle  between  helicopter  velocity  vector  VI NIT  and  local 
horizontal 

DNSCH  -  Dynamic  pressure  scaler 

DP  IT  -  Offset  in  pitch  of  motor  nozzle 

DT  -  Delta  time  for  integration  routine  after  TDC 

DYAW  -  Offset  in  yaw  of  motor  nozzle 

DZZ  -  Z  component  of  downwash  force 

EOL  -  End  of  launch  flag 

FACTOR  -  Variable  associated  with  motor  nozzle  (ft.) 

FCA  -  Forward  cant  angle  of  helicopter 
FTORM  -  Output  format  flag  (ft.  or  meters) 

GREFF  -  Force  of  gravity  (ft. /sec*) 

HELIC  -  Flag  for  use  of  downwash  profile  (l.-use,0.-do  not  use) 

PITMAL  -  Pitch  mal launch  rate 

QELNCH  -  Uauncher  elevation  angle  with  respect  to  helicopter  centerline 
REF  -  Radius  of  the  earth 

RK  -  Flag  for  type  of  Runge-Kutta  integration  routine  (2-2nd  order ,4-4th 
order). 

ROTRAD  -  Radius  of  helicopter  rotor  blades 

SCARF  -  Variable  associated  with  motor  nozzle  (ft.) 

TDC  -  Time  that  integration  delta  time  DELTT  changes  to  DT. 

TFINAL  -  Maximum  time  of  flight 

TIME  -  Program  counter 

TLL  -  Time  missile  left  launcher 

TPRINT  -  Time  increment  for  printing  results 

TROT  -  Time  of  rotation  (sec.) 

VINIT  -  Intial  helicopter  velocity  vector's  magnitude  (airspeed  ft. /sec. 
WCF  -  Crossrange  wind  force 
WDF  -  Downrange  wind  force 

XMV-X  -  Component  of  missile  velocity  at  launch  (ft. /sec.) 

XTAR  -  Range  of  target  (ft.  Earth  fixed  coordinate  system) 

YAWMAL  -  Yaw  mal launch  rate 

YTAR  -  Crossrange  of  target  (ft.  Earth  fixed  coordinate  system) 

ZTAR  -  Altitude  of  target  (ft.  Earth  fixed  coordinate  system) 

INTERGERS 

N-  Integration  routine  upper  limit 

NM  -  Integration  routine  lower  limit _ 
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2.  Computed  Variables 


Table  X  is  a  list  of  variables  computed  in  the  2.75  digital 

simulation. 

TABLE  X.  2.75  COMPUTED  VARIABLES 


Logical? 

BURNOUT  -  Motor  has  burned  out 

EOL  -  [END  OF  LAUNCH]  Trajectory  begins  with  parameters  at  time  rocket 
clears  tube 

IMPACT  -  Impact  occurred 

LNEXIT  -  Rocket  clear  of  launcher 

LRATE  -  Rocket  CG  has  cleared  launcher 

NOTIME  -  Runtime  exceeds  allowed  TFINAL 

SETBO  -  Set  rocket  parameters  to  burnout  values 

SKIP  -  Overrides  SETBO 

Real  s 

1X0  -  Initial  axial  moment  of  inertia  (slugs/ft2) 

IXF  -  Final  axial  moment  of  inertia  (slugs/ft2) 

IYO  -  Initial  transverse  moment  of  inertia  (slugs/ft2) 

IYF  -  Final  transverse  moment  of  inertia  (slugs/ft2) 

IX  -  Axial  moment  of  inertia 
IY  -  Transverse  moment  of  inertia 
IZ  -  Vertical  moment  of  inertia 

IXDOT  -  Time  rate  of  change  of  axial  moment  of  inertia  (slugs. sec/ft2) 

IYDOT  -  Time  rate  of  change  of  transverse  moment  of  inertia  (slugs. sec/ft2) 

IZDOT  -  Time  rate  of  change  of  vertical  moment  of  inertia  (slugs. sec/ft2) 

ISP  -  Specific  impulse 

MCDOFF  (21)  -  Mach  number  table  for  coefficients  of  DRAG  w/power  off 

MCDON  (21)  -  Mach  number  table  for  coefficients  of  DRAG  w/power  on 

MCMQ  (18)  -  Mach  number  table  for  pitch  damping  coefficients 

MCNA  (7)  -  Mach  number  table  for  normal  force 

MCP  (9)  -  Mach  number  table  for  center  of  pressure 

MOOT  -  Mass  flow  rate  as  a  function  of  thrust 

ACTI  (3,3)  -  Inverse  of  DCM  Matrix 

CDPOFF  (21)  -  Coefficients  of  DRAG  with  power  off  (function  of  Mach  data) 
CDPON  (21)  -  Coefficients  of  DRAG  with  power  on  (function  of  Mach  data) 

CMQT  (17)  -  Coefficients  of  pitch  damping  (rad/sec  -  function  of  mach  data) 
CNAT  (7)  -  Slope  coefficient  of  normal  force  (f/deg.  -  function  of  mach 
data) 

CPT  (8,2)  -  Center  of  pressure  (calibers  from  nose  -  function  of  mach  data 
and  angle  of  attack) 

CPAF  (2)  -  Angle  of  attack  table  (degrees) 

CK  (4,  50)  Calculated  integrating  constants 
LABEL  (8)  -  Input  of  trajectory  characteristics 
IDAY  (3)  -  Date  (internal  to  computer) 

DCM  (3,3)  -  Euler  angle  transformation  matrix 
DTDH  (12)  -  Atmospheric  density  reference  table 
HGEOPB  (13)  -  Atmospheric  altitude  reference  table 
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TABLE  X.  2.75  COMPUTED  VARIABLES  (Cont Jd) 

*  — - -  — "  ■ 

;  PRESB  (12)  -  Atmospheric  pressure  reference  table 

SPINT  (30)  -  Time  breakpoints  for  SPNRAT  (table) 

|  SPNRAT  (30)  -  Average  roll  rate  (rev/sec  -  table) 

TEMPMB  (12)  -  Atmospheric  temperature  reference  tables 
j  THRUST  (9)  -  Thrust  table  (pounds) 

»  TTIME  (9)  -  Time  table  for  thrust 

J  WUVW  (3,3)  -  Result  of  DCM  (3,3)  *  WXYZ  (3,3)  (Effect  of  wind  on  rocket 

1  body) 

WXYZ  (3,3)  -  Coefficients  of  wind 

XX  (3,50)  -  XX  (1,50)  Derivatives  of  rocket  trajectory  parameters 
XX  (2,50)  State  of  rocket  trajectory  parameters 
XX  (3,50)  Not  used 

YY  (3,50)  -  YY  (1,50)  Derivative  at  beginning  of  time  step 
YY  (2,50)  State  at  beginninq  of  time  step 
YY  (3,50)  Not  used 
CLP  -  Roll  damping  coefficient 

!  CG  -  Instantaneous  center  of  gravity  (ft.  from  tail) 

!  CG0  -  Initial  center  of  gravity  (ft.  from  tail) 

CGF  -  Final  center  of  gravity  (ft.  from  tail) 

DREF  -  Reference  diameter  of  rocket  (ft) 

F  -  Air  density  constant 
PI  -  Variable  equivalent  to  PI  =  3.14159 
|  RADCON  -  Degrees  per  radian  =  57.2957795 

!  S  -  Surface  area  of  earth 

SUSOFF  -  Motor  burn  time 
TI  -  Total  impulse 
WO  -  Initial  rocket  weight 
WF  -  Final  rocket  weight 
WL  -  Weight  when  rocket  clears  launcher 
XLENG  -  Missile  length 

CGDOT  -  Time  rate  of  change  of  center  of  gravity 
XMASS  -  Initial  mass  of  rocket 
QEMILS  -  Quadrant  elevation  (mils) 

AIRSP  -  Airspeed  (knots) 

HELAT  -  Helicopter  attitude  (degrees) 

AREFF  -  Reference  area  of  rocket  (ft2) 

THL  -  Angle  between  launcher  centerline  and  local  horizontal 

ALV  -  Angle  between  launcher  centerline  and  helicopter  velocity  vector  VINIT 

TTO  -  Time  till  thrust  equals  zero 

TEX  IT  -  Time  missile  exits  launcher 

DXLD  -  X  -  Coordinate  of  launcher  CG  with  respect  to  the  helicopter  rotor 
hub 

DYLN  -  Y  -  Coordinate  of  launcher  CG  with  respect  to  the  helicopter  rotor 
hub 

DZLN  -  Z  -  Coordinate  of  launcher  CG  with  respect  to  the  helicopter  rotor 
hub 

DXX  -  X  -  Components  of  downwash  coefficient 

DXY  -  Y  -  Components  of  downwash  coefficient 

DZZ  -  Z  -  Components  of  downwash  coefficient 

SI  -  Euler  Angle  PSI 
TH  -  Euler  Angle  THETA 

FI  -  Euler  Angle  PHI _ 
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TABLE  X.  2.75  COMPUTED  VARIABLES  (Cont'd) 


SSI  -  SIN  (PSI) 

CSI  -  COS  (PSI) 

STH  -  SIN  (THETA) 

CTH  -  COS  (THETA) 

SFI  -  SIN  (PHI) 

CFI  -  COS  (PHI) 

P  -  Angular  acceleration  about  rocket  Z-axis 

UDOT  -  Linear  acceleration  along  rocket  X-axis 

VDOT  -  Linear  acceleration  along  rocket  Y-axis 

WDOT  -  Linear  acceleration  along  rocket  Z-axis 

CMQ  -  Rate  Damping  coefficient  for  rates  about  rocket 

AUX  -  Variable  equal  to  (Dynamic  Pressure  *  rocket  diameter  squared) 

2  *  Relative  Velocity 

XMP1AS  -  Aerodynamic  torques  about  rocket  X-axis 

XMQ1AS  -  Aerodynamic  torques  about  rocket  Y-axis 

XMR1AS  -  Aerodynamic  torques  about  rocket  Z-axis 

CP  -  Center  of  pressure  as  a  function  of  Mach  No.  and  Angle  of  Attack 
AMA  -  Aerodynamic  moment  arm 
SM  -  Stability  margin  (calibers) 

XMAY  -  Aerodynamic  restoring  torques  about  rocket  Y-axis 

XMAZ  -  Aerodynamic  restoring  torques  about  rocket  Z-axis 

TAX  -  Torque  about  rocket  X-axis  due  to  angular  velocities 

TAY  -  Torque  about  rocket  Y-axis  due  to  angular  velocities 

TAZ  -  Torque  about  rocket  Z-axis  due  to  angular  velocities 

TIX  -  Torque  about  rocket  X-axis  due  to  time  rate  of  change  of  inertia 

TIY  -  Torque  about  rocket  Y-axis  due  to  time  rate  of  change  of  inertia 

TIZ  -  Torque  about  rocket  Z-axis  due  to  time  rate  of  change  of  inertia 

XMPIGS  -  Torque  about  rocket  X-axis  due  to  components  of  thrust 

ARM  -  Length  of  moment  arm 

XMQ1GS  -  Torque  about  rocket  Y-axis  due  to  components  of  thrust 

XMR1GS  -  Torque  about  rocket  Z-axis  due  to  components  of  thrust 

TORQX  -  Total  torque  about  rocket  X-axis 

TORQY  -  Total  torque  about  rocket  Y-axis 

TORQZ  -  Total  torque  about  rocket  Z-axis 

ROLL  -  Roll  rate  (hz) 

GAX  -  Gravitational  acceleration  along  rocket  X-axis 

GAY  -  Gravitational  acceleration  along  rocket  Y-axis 

GAZ  -  Gravitational  acceleration  along  rocket  Z-axis 

UDOTL  -  Acceleration  rocket  along  X-AXIS  due  to  drag,  thrust  +  gravity  only 

UL  -  Velocity  of  rocket  along  launcher  centerline 

REZIF  -  Rocket  altitude  from  center  of  earth 

RF  -  Distance  from  center  of  earth  to  launcher 

HEIGHT  -  Distance  from  earth's  crust  to  launcher 

SQIRT  -  Horizontal  distance  from  launcher  to  rocket 

DYNAR  -  Dynamic  pressure  times  reference  area  (ft) 

DZRFF  -  Reference  diameter  (ft) 

TANG  -  Cosine  of  trajectory  angle  from  center  of  earth 
RGF  -  Gravitation  reference  height  (ft) 

WIND  -  Magnitude  of  wind  force 
TPHASE  -  (TIME  +  TROT  =  .187) 

XH  -  Instantaneous  distance  from  rotor  hub  to  rocket  along  X-axis _ 


18 


TABLE  X.  2.75  COMPUTED  VARIABLES  (Cont 'd) 


YH  -  Instantaneous  distance  from  rotor  hub  to  rocket  along  Y-axis 
ZH  -  Instantaneous  distance  from  rotor  hub  to  rocket  along  Z-axis 
XMXH  -  Initial  position  of  rocket  with  respect  to  rotor  hub  along  X-axis 
YMYH  -  Initial  position  of  rocket  with  respect  to  rotor  hub  along  Y-axis 
ZMZH  -  Initial  position  of  rocket  with  respect  to  rotor  hub  along  Z-axis 
DDIS  -  Vertical  position  of  rocket  with  respect  to  launcher 
THAT  -  Total  Horizontal  Attitude  of  Helicopter  (deg) 

THEO  -  Downrange  rocket  trajectory  with  respect  to  vertical  (deg) 

THTH  -  Downrange  rocket  trajectory  with  respect  to  horizontal  (deg) 

XD  -  Horizontal  position  of  rocket  with  respect  to  rotor  hub 

YD  -  Deflection  position  of  rocket  with  respect  to  rotor  hub 

ZD  -  Vertical  position  of  rocket  with  respect  to  rotor  hub 

DDISS  -  Horizontal  distance  from  rocket  to  rotor  hub 

WXYZ  (1,1)  -  Magnitude  of  wind  X-axis  vector 

WXYZ  (2,1)  -  Magnitude  of  wind  Y-axis  vector 

WXYZ  (3,1)  -  Magnitude  of  wind  Z-axis  vector 

HGTIND  -  Height  indicator 

HGEOP  -  Atmospheric  altitude 

WMOL  -  Molecular  weight 

TEMPMD  -  Instantaneous  observed  atmospheric  temperature 
HGEOPO  -  Instantaneous  observed  atmospheric  altitude 
DODH  -  Instantaneous  observed  atmospheric  density 
PRESO  -  Instantaneous  observed  atmospheric  pressure 
TEMPM  -  Median  temperature  (molecular  temperature) 

TEMP  -  Actual  temperature 

VSOUND  -  Velocity  of  sound  (function  of  altitude,  temperature,  &  air 
denstiy) 

PRES  -  Air  pressure 

XMW  -  Apparent  velocity  of  rocket  in  X  vector  (m-wind) 

YMW  -  Apparent  velocity  of  rocket  in  Y  vector 
ZMW  -  Apparent  velocity  of  rocket  in  Z  vector 
VVW  -  Vertical  velocity  (squared) 

VRF  -  Relative  velocity  of  rocket  (ft/s) 

XMACHN  -  Relative  velocity  (MACHS) 

VALF  -  Slope  of  velocity  vector 
WUVW  (1,1)  -  Wind  X  vector  on  rocket 
WUVW  (2,1)  -  Wind  Y  vector  on  rocket 
WUVW  (3,1)  -  Wind  Z  vector  on  rocket 
ALFTOF  -  Total  angle  of  flight 

BANK  -  Angle  of  heliocopter  body  with  respect  to  vertical 
ALPPRN  -  Vertical  body  angle  with  respect  to  downrange 
BETPRN  -  Horizontal  body  angle  with  respect  to  downrange 
CN  -  Normal  force  coefficient 
CY  -  Side  force  coefficient 
CA  -  Axial  force  coefficient 
THRSTP  -  Magnitude  of  thrust  vector 
DP  -  Dynamic  pressure 

DPS  -  Dynamic  pressure  on  surface  of  missile 

FU1AP  -  Aerodynamic  force  along  rocket  X-axis 

FV1AP  -  Aerodynamic  force  along  rocket  Y-axis 

FW1AP  -  Aerodynamic  force  along  rocket  Z-axis _ 


TABLE  X.  2.75  COMPUTED  VARIABLES  (Cont'd) 


FU1GP  -  Thrust  force  along  rocket  X-axis 

FV1GP  -  Thrust  force  along  rocket  Y-axis 

FW1GP  -  Thrust  force  along  rocket  Z-axis 

UMF1  -  Linear  acceleration  along  rocket  X-axis 

VMF1  -  Linear  acceleration  along  rocket  Y-axis 

WMF1  -  Linear  acceleration  along  rocket  Z-axis 

U  -  Linear  velocity  along  rocket  X-axis 

V  -  Linear  velocity  along  rocket  Y-axis 

W  -  Linear  velocity  along  rocket  Z-axis 

P  -  Angular  velocity  about  rocket  X-axis 

Q  -  Angular  velocity  about  rocket  Y-axis 

R  -  Angular  velocity  about  rocket  Z-axis 

PDOT  -  Angular  acceleration  about  X-axis 

QDOT  -  Angular  acceleration  about  Y-axis 

RDOT  -  Angular  acceleration  about  Z-axis 

PNOM  -  Length  of  rocket  vertical  velocity  vector 

TWOPI  -  Two  times  the  radian  angle  PI  (Circumference  of  unit  circle) 

PHI  I  -  Angle  PHI 

SNPHI  -  Sine  of  PHI  I 

CMP HI  -  Cosine  of  PHI  I 

ETAPRN  =  Angle  between  horizontal  and  bank  vector 

ALPHPR  -  Angle  Alpha  (Rad)  -  Angle  between  Z  and  X  velocity  vectors 

BETAPR  -  Angle  Beta  (Rad)  -  Angle  bewteen  Y  and  X  velocity  vectors 

SIDOT  -  Derivative  of  PSI 

THDOT  -  Derivative  of  THETA 

FIDOT  -  Derivative  of  PHI 

TPR  -  Time  till  print  (Run  time  [TIME  -  TPRINT]  Time  of  period) 

ALPTD  -  Angle  Alpha  in  degrees 
BETTD  -  Angle  Beta  in  degrees 
ETATD  -  Angle  Eta  in  degrees 
ALF  -  Angle  TTIpTia  in  radians 
BETO  -  Angle  Beta  in  radians 
TP1  -  Time  (counter  1) 

XP1  -  Rocket  downrange  position  (time  1) 

YP1  -  Rocket  cross-range  position  (time  1) 

ZP1  -  Rocket  altitude  position  (time  1) 

XP2  -  Rocket  downrange  impact 
YP2  -  Rocket  cross-range  impact 

ZP2  -  Rocket  altitude  (impact  -  ground  coordinate  system) 

TP2  -  Time  rocket  impact 

XMP  -  Impact  range 

TMP  -  Time  of  impact 

PID  -  Roll  rate  (rev/sec)  at  impact 

WINX  -  Velocity  along  X-axis  including  wind 

including  wind 


WINY  -  Velocity  along  Y-axis 


-  Velocity  along  Z-axis  including  wind 
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TABLE  X.  2.75  COMPUTED  VARIABLES  (Cont'd) 


ALFO  -  Vertical  angle  of  impact 
QID  -  Pitch  at  impact  (Radians) 

RID  -  YAW  at  impact  (Radians) 

DIST  -  Distance  to  go  (rocket  to  target) 

PID1  -  Angular  acceleration  about  rocket  X-axis  (degrees) 

Q1D1  -  Angular  acceleration  about  rocket  Y-axis  (degrees) 

R1D1  -  Angular  acceleration  about  rocket  Z-axis  (degrees) 

KN  -  Integration  counter 
CK  (1,4)  -  Integration  coefficient  array 


3.  Outputs 
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There  are  two  output  sources  available  to  the  user  of  the  2.75 
digital  simulation.  The  first  is  output  by  the  program  via  the  line  printer. 
The  second  is  output  to  a  local  file  for  use  with  the  plotting  routine  associ¬ 
ated  with  this  program. 

(b)  Printed  Variables 

There  are  three  sets  of  outputs  to  the  line  printer.  The 
first  is  the  set  of  all  initial  conditions  for  the  simulation.  The  second  set 
is  output  after  the  value  of  counter  exceeds  the  value  of  a  print  flag.  The 
last  set  is  output  after  impact  or  when  the  simulation  runs  past  the  maximum 
allowable  time.  A  listing  of  these  variables  and  their  definitions,  in  the 
order  they  are  output,  follows  in  Table  XI. 


TABLE  XI.  2.75  PRINTED  OUTPUT  VARIABLES 


T I  ME!  -  Time  of  Flight 
RANGE  -  Downrarige 
DEFL  -  Crossrange 
ALT  -  Altitude 

VELR  -  Relative  Velocity  (Magnitude) 

THRUST  -  Magnitude  of  Thrust 

XDOT  -  Linear  Velocity  Along  Ground  X{ Range) -Axis 
YDOT  -  Linear  Velocity  Along  Ground  Y(Deflection)-Axis 
ZDOT  -  Linear  Velocity  Along  Ground  Z(Vertical  )-Axis 
MACH  -  Mach  Speed 
MASS  -  Mass  of  Rocket 

UUU  -  Total  Linear  Velocity  Along  Rocket  X-Axis 

VVV  -  Total  Linear  Velocity  Along  Rocket  Y-Axis 

WWW  -  Total  Linear  Velocity  Along  Rocket  Z-Axis 

MOOT  -  Time  Rate  of  Change  of  Rocket  Mass 

(SM)  - 

UDOT  -  Total  Linear  Acceleration  Along  Rocket  X-Axis 
VDOT  -  Total  Linear  Acceleration  Along  Rocket  Y-Axis 
WDOT  -  Total  Linear  Acceleration  Along  Rocket  Z-Axis 
DP  -  Dynamic  Pressure 

RHO  -  Air  Density _ 
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TABLE  XI.  2.75  PRINTED  OUTPUT  VARIABLES  (Cont'd) 


P  -  {Roll  Rate)  Angular  Velocity  About  Rocket  X-Axis 
Q  -  (Pitch  Rate)  Angular  Velocity  About  Rocket  Y-Axis 
R  -  (Yaw  Rate)  Angular  Velocity  About  Rocket  Z-Axis 
CA  -  Axial  Force  Coefficient 
ALFTOT  -  Total  Angle  of  Flight 

PID1  -  Angular  Acceleration  About  Rocket  X-Axis  (Degrees) 

QID1  -  Angular  Acceleration  About  Rocket  Y-Axis  (Degrees) 

RID1  -  Angular  Acceleration  About  Rocket  Z-Axis  (Degrees) 

CN  -  Vertical  Coefficient  of  Force 
ALPTD  -  Angle  Alpha  (Degrees) 

TORQX  -  Total  Torque  ^bout  Rocket  X-Axis 
TORQY  -  Total  Torque  About  Rocket  Y-Axis 
TORQZ  -  Total  Torque  About  Rocket  Z-Axis 
CY  -  Side  Force  Coefficient 
BETTD  -  Angle  Beta  (Degrees) 

WINX  -  Velocity  ALong  X-Axis  Including  Effect  of  Wind 

WINY  -  Velocity  Along  Y-Axis  Including  Effect  of  Wind 

WINZ  -  Velocity  Along  Z-Axis  Including  Effect  of  Wind 

LNEXIT  -  Logical  Representing  Rocket  Clear  Launcher 

XMQ1AS  -  Aerodynamic  Torque  About  Rocket  Y-Axis 

XMAY  -  Aerodynamic  Restoring  Torque  About  Rocket  Y-Axis 

TAY  -  Torque  About  Rocket  Y-Axis  due  to  Angular  Velocities 

TIY  -  Torque  About  Rocket  Y-Axis  due  to  Time  Rate  of  Change  of  Inertia 

CMQ  -  Rate  Damping  Coefficient  for  Rates  About  Rocket  Axes 

XMR1AS  -  Aerodynamic  Torque  About  Rocket  Z-Axis 

XMAZ  -  Aerodynamic  Restoring  Torque  About  Rocket  Z-Axis 

TAZ  -  Torque  About  Rocket  Z-Axis  due  to  Angular  velocities 

TIZ  -  Torque  About  Rocket  Z-Axis  due  to  Time  Rate  of  Change  of  Inertia 

BANK  -  Angle  Between  Positive  Vertical  and  Helicopter  Vertical  Centering 

PSI  -  Euler  Angle  psi 

THETA  -  Euler  Angle  theta 

PHI  -  Euler  Angle  phi 

UMF1  -  Linear  Acceleration  Along  Rocket  X-Axis 
VMF1  -  Linear  Acceleration  Along  Rocket  Y-Axis 
WMF1  -  Linear  Acceleration  Along  Rocket  Z-Axis 
GAX  -  Gravitational  Acceleration  Along  Rocket  X-Axis 
GAY  -  Gravitational  Acceleration  Along  Rocket  Y-Axis 
GAZ  -  Gravitational  Acceleration  Along  Rocket  Z-Axis 

CGTRVL  -  Distance  Center  of  Gravity  of  Rocket  has  Travelled _ 

(b)  Plot  Tape  Variables 

The  following  is  a  list  of  the  2.75  variables  that  are 
written  to  the  plot  tape  from  the  6-DOF  simulation.  These  are  reproduced  in 
the  order  that  they  are  stored  on  record. 
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TABLE  XII.  2.75  PLOT  TAPE  VARIABLES 


Time  of  FI i ght  ( Sec . ) 

Range 

Deflection  -  Positive  Right 
Altitude  -  Positive  Down 
Relative  Velocity 
Mass  (Slugs) 

Range  Velocity 

Cross  Range  Velocity  -  Postivie  Right 
Vertical  Velocity  -  Positive  Down 
Mach  Number 
Thrust 

Roll  Rate  (Rev/Sec) 

Pitch  Rate  (Deg/Sec) 

Yaw  Rate  (Deg/Sec) 

Dynamic  Pressure 
Range  to  go  to  Target 
Roll  Acceleration  (Rad/Sec^l 
Pitch  Acceleration  (Rad/Sec^) 

Yaw  Acceleration  (Rad/Sec^) 

(HALFI) 

Center  of  Gravity  -  Ft.  From  Tail 
Rocket  X-Axis  Velocity  -  U 
Rocket  Y-Axis  Velocity  -  V 
Rocket  Z-Axis  Velocity  -  W 
Densi  ty 

Center  of  Pressure  -  Ft.  From  Tail 
Body  X-Axis  Acceleration  -  UDOT 
Body  Y-Axis  Acceleration  -  VDOT 
Body  Z-Axis  Acceleration  -  WDOT 
(ALFTO) 

(PMOMA) 

Angle  ALPHA  (Degrees) 

Angle  BETA  (DEG) 

Angle  ETA  (DEG) 

(FORSU) 

(QMONA) 

Angle  THETA  (DEG) 

Angle  PHI  (DEG) 

Angle  PSI  (DEG) 

(FORSV) 

(RMONA) 

Relative  Range  Velocity 

Relative  Deflection  Velocity  -  Positive  Right 
Relative  Vertical  Velocity  -  Positive  Down 
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TABLE  XII.  2.75  PLOT  TAPE  VARIABLES  (Cont’d) 


TFO'R'SW) 

( PMOMG ) 

Wind  Force  Along  Range  Axis 
Wind  Force  Along  Deflection  Axis 
Wind  Force  Along  Vertical  Axis 
Time  Rate  of  Change  in  Mass  (Slugs/Sec) 

( QMOMG )  - 
{ RMOMG )  - 

Axial  Coefficient  of  Force 
Side  Coefficient  of  Force 
Vertical  Coefficient  of  Force 
(NALP) 

{ BETI ) 

Bank  Angle 
( HCL ) 

Rate  Damping  Coefficient  for  Rates  About  Rocket  Axes. 


V.  TEKPLOT  PLOTTING  PROGRAM 

A.  Introduction 

The  2.75  digital  simulation  takes  advantage  of  a  multipurpose  plot 
routine  developed  at  MICOM.  This  program,  TEKPLOT,  enables  the  user  to  quickly 
generate  readable  and  understandable  evaluation  packages  for  each  data  run. 

In  association  with  TEKPLOT  is  a  routine,  RPLOT ,  which  ensures  that  the  data 
input  to  the  TEKPLOT  routine  is  in  a  suitable  format. 

B.  RPLOT  Tape  Read/Write  Program 

RPLOT  is  a  program  designed  to  write  data  from  an  unformatted  plot 
file  to  a  formatted  tape.  Its  listing  is  presented  in  Appendix  E. 

C.  TEKPLOT  Program  Operations 

TEKPLOT  is  a  program  designed  to  generate  x-y  plots  at  a  TEKTRONIX 
4014  terminal,  where  hardcopies  can  be  made.  Its  listing  is  presented  in 
Appendix  F.  The  program  is  designed  to  run  interactively,  using  a  local  file 
as  the  data  file.  Normally,  TAPE2  is  used  as  the  local  data  file  when  the 
following  conditions  are  met. 

1.  The  number  of  data  points  for  each  variable,  NPTS,  is  on  the 
first  record  of  TAPE2. 

2.  The  minimum  and  maximum  values  for  each  variable  are  stored  in 
the  next  NV  records  where  NV  is  the  number  of  variables. 

These  conditions  are  present  when  the  output  from  2.75  is  in  its  standard 
form.  If  both  conditions  are  not  met,  the  format  is  non-standard  and  the 
operator  must  answer  a  displayed  question  regarding  the  conditions  with  a 
"NO".  At  this  time,  the  program  reads  through  TAPE2  and  finds  the  number  of 
data  points  for  each  variable  and  the  minimum  and  maximum  values  of  each 


variable.  These  values  are  written  to  TAPE10  in  standard  format.  NPTS  blocks 
of  records  are  then  written  to  TAPE10  with  each  record  holding  the  first, 
second,  etc.  data  point  of  all  records.  This  process  of  recording  is  done 
by  SUBROUTINE  FORMAT. 

This  program  assumes  that  there  are  NV  variables  on  each  record  of  the 
data  tape  and  that  the  operator  knows  the  sequential  order  of  the  variables. 

The  operator  has  the  option  of  choosing  any  pair  of  variables  to  be  plotted 
and  can  designate  either  of  the  variables  to  be  the  independent  variable.  The 
program  provides  the  option  for  plotting  1  to  N  pairs  of  variables  at  a  time, 
with  a  maximum  N  of  60. 

Under  current  system  restraints  the  program  enters  a  pause  mode  after 
plotting  to  allow  for  a  hard  copy  of  the  screen.  To  proceed  to  the  next  plot 
the  operator  must  enter  any  character  from  the  keyboard.  Sometimes  it  may  be 
necessary  to  send  the  cursur  home  and  then  enter  any  character  from  the  key¬ 
board.  The  program  assumes  that  the  current  data  tape  contains  one  run  of 
data.  At  the  end  of  each  series  of  plots  the  operator  is  given  the  option  of 
whether  or  not  he  wants  to  continue  plotting  from  the  current  data  tape.  If 
a  "NO"  is  entered,  the  program  stops  execution.  If  "NO"  is  not  entered,  the 
operator  can  re-initialize  the  program  and  continue  to  plot  from  the  current 
data  tape.  If  he  chooses  to  plot  from  the  current  data  tape,  he  will  use 
TAPE2  if  the  original  data  was  in  standard  format.  If  the  data  was  in  non¬ 
standard  format,  he  should  return  TAPE2  and  rename  TAPE10  as  TAPE2.  If  the 
operator  wants  to  plot  from  a  different  data  run,  TAPE2  and  TAPE10,  if 
applicable,  must  be  returned.  Then  he  must  attach  or  create  a  new  TAPE2  and 
proceed  to  re-initialize  the  program.  If  one  had  a  file  called  DATA  that 
contained  several  runs,  each  followed  by  an  END  OF  RECORD,  the  tenth  run  could 
be  plotted  by  copying  records  onto  a  file  called  DUMMY.  The  tenth  record  could 
then  be  copied  into  a  file  called  TAPE2.  TAPE2  rewound,  and  the  program 
initialized. 

SUBROUTINE  LABLR  creates  tables  for  the  x  and  y  axis  from  a  block  data 
statement  called  LAB.  The  program  is  designed  to  handle  60  labels  with  each 
label  containing  a  maximum  of  3  words.  The  labes.,  in  block  data,  appear  in 
the  same  sequential  order  that  the  variables  appear  on  the  data  tape.  If  the 
operator  wants  to  change  the  label  of  the  Ith  variable  on  the  record,  he  simply 
chang~s  the  Ith  label  in  block  data.  The  program  will  correctly  select  the 
labels  for  the  specified  pair  of  variables  to  be  plotted. 

SUBROUTINE  XYGRID  draws  the  grid  for  the  XY  plot.  If  desired,  you  can 
change  the  tic  mark  form  of  the  grid  by  changing  this  subroutine. 

SUBROUTINE  FORMAT  changes  data  received  in  non-standard  form  to  standard 
form  (described  above) . 

VI.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  study  of  PERSHING  II  impact  points  has  been  successfully  completed  for 
important  debris  through  the  use  of  the  six  degree  of  simulation  program 
MSIXB1.  The  program  has  been  described  so  that  it  may  be  used  as  a  tool  to 
aid  in  future  studies.  Important  subroutines  including  ballistic  and  atmos¬ 
pheric  coding  have  been  added  and  successfully  employed  in  MSIXB1.  The  2.75 
simulation,  which  may  also  be  used  for  debris  studies,  has  been  described. 


Thus,  the  PERSHING  II  debris  dynamics  simulation  study  has  been  successful, 
not  only  in  providing  needed  data,  but  also  in  providing  tools  for  the  future. 

Recommendations  for  the  MSIXB1  model  are  to  further  investigate  the 
integration  routine  and  the  quaternion  transformations.  When  high  spin  rate 
debris  from  PERSHING  II  was  investigated,  the  integration  result  varied  widely 
as  the  integration  step  A(30)  was  increased.  Extensive  testing  of  the  program 
under  high  spin  rate  conditions  using  various  integration  steps  was  performed, 
and  a  .0008  second  step  was  determined  to  be  the  maximum  valid  step  size.  Any 
larger  size  caused  a  response  in  which  the  impact  piont  became  continuously 
further  from  the  precise  point  calculated  using  smaller  step  sizes.  Also,  a 
complete  documentation  of  how  MSIXBl's  quaternion  equations  are  developed  from 
general  quaternion  equations  would  be  helpful  to  users. 

Of  interest  to  users  who  use  multiple  runs,  changing  inputs  under  program 
control  each  time  would  be  a  statistical  extension  of  the  program.  This  could 
be  accomplished  by  writing  values  to  a  data  file  as  they  are  calculated  by  the 
program.  This  file  could  then  be  accessed  by  a  statistical  routine.  This 
process  is  now  being  developed  to  analyze  impact  data  using  a  circular  error 
pattern  statistical  routine. 
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SBATCH 
SPROG  MAIN 

PROGRAM  MSIXB1 

C  *********  6DOF  SIMULATION  ************* 

C  THIS  PROGRAM  WAS  ORIGINALLY  WRITTEN  BY  MARTIN  MARIETTA  C 

C  IT  HAS  BEEN  MODIFIED  BY  PATRICK  TI L LEY . GERAL DN JOHNSON ,  C 

C  DONN  HALL  AND  OTHERS  AT  UAH  TO  ACCOMODATE  DEBRIS  STUDIES  C 

LOGICAL  PFLAG.QRFLAG 
REAL  Ml, M2, M3 
EXTERNAL  SIGN 
COMMON/DATA/  A{  2 8j0Tj0T > 

COMMON  /ALLT/  DELT , TIME , IT IME , DELT2 

COMMON  /AKIN/ADR.ADP ,ADY,ADF1 , ADF 2 , ADF 3 , ADM  1 , ADM2 , ADM3 ,AWX ,AWY , AW2 

COMMON  /MONT /NMRC , LRC , NRC , IB 

COMMON/KILL/  PE . CEP50 , THETAA , S 1 GTHE , A249 , PK 

COMMON  /AKOT/AA1 , AA2 , AA3 . AAX . AAY , AA2 , AVX . AVY , AV2 , AX , AY , AZ , AP D , 

1  AQD , ARD , AP , AQ , AR , APS  I , ATHETA. APH I , AVA , AVT .AMACH , ALPHA, AB ETA ,AETA, 

2  AQDP,  AH.  AW , ACG , A IXX , AI YY , ACA , ACMT , ACL  P , ACL  D , ACL  1 , 

3  ACMP , ACMY , ACNP , ACNY , AC  1 X , AC2X , AC3X , AC  1 Y , AC2 Y , AC3Y , AC  1 Z , AC2Z , AC3Z 
4 , APH 1 D , ATHED .APS  ID 

COMMON /GNC/ACQ , 1 E ,  ESP , ESY . GP , GY , PRS , DBP , DBY , XP , YP , ZP , 

1  ACQRG.ACQALT.TSEQ,  SEARCH , TPR 1  NT , P HASE , CPS , CYS , TSRCH , GGU I D 

COMMON  /PULSES/  NPLS!360) 

COMMON  /TARG/XBGO, YBGO ,XFGO , YFGO , NBG , NFG , NTRUE , NTOT , OPLS 
COMMON  /CAERO/  CNA , CND , CMA , CMD . V23S , V23 , CAFAC , 

1  AVI , AV2 , AV3 , AV4 , ABSA,ABSB,SAD,SBD,XX,CL, 

2  DPS,DYS,DCY,DCG,D2V,V2A,V3A,SA.SB,F1,F2,F3,M1,M2,M3 
COMMON  /SPLOT/  QRMX , D I VE , TURN , DEL  TAX , DEL TAY 

COMMON  /CPLOT/  HERR , DFL AR E , HGO , ADPX , AD YX , EL  OS 
COMMON  /SAUTO/  QFB , RFB , DTHS , R 
DIMENSION  SYM ( 27  > 

CHARACTER*70  TITLE 
CHARACTER  *3  STA 
CHARACTER  *12  F NAME 
CHARACTER  *9  FM 

DATA  SYM/’TIME' , 'PSI  ',’PH1  ’.’THET’.’P  ’,’Q  ’,’R 

*  ’PDOT’ , *QDOT’ , 'RDOT' , 'AX  ’ , ’AY  ’ , ’AZ  ’,’M1  ’, 

*  ’M2  ’ , 'M3  ' , ’ XM  ’ , ’ YM  ’,’ZM  ’,’VXM  ’ , ’ VYM  ’, 

*  ’  VZM  ’ , ’DLAX’ , ’DLAY’ , 'FI  ’,’F2  ’,’F3  ’/ 

OPEN! 10,FILE-’RUN.DAT’ , STATUS- ’OL D ’ ) 

READ! 10,807)  NRUN 
807  FORMAT! 15) 

NRUN-NRUN+ 1 
CLOSE! 10) 

OPEN! 10, FILE- ’RUN. DAT’ , STATUS- ’ OLD ’ > 

WRITE! 10, B07)  NRUN 

DATA  DPR, ZERO, ZN I NE/57, 2957795,0. ,999./ 

515  IB-0 

COMMENT  -  INITIAL  CALCULATIONS 
1  CALL  I NP  UT { A , 2B00, 300, KRN ) 

IF! IB.EQ. 1 >  GOTO  515 
IF!  NMRC.LE.0  )  WRITE  !3,100)  KRN 
IF!  NMRC.LE.0  )  WRITE  16,100)  NRUN 
100  FORMAT  { IX, ’RUN’ , 15) 

IF!  NMRC.GT.0  >  WRITE  13,1100)  KRN, NMRC 
IF!  NMRC.GT.0  )  WRITE  16.1100)  KRN , NMRC 
1100  FORMAT !4H1RUN,I5,2X, 2HMC , 15/ > 

WR1TE16, 1001 > 

IF!  A! 25 ) . GT .0 .  )  WRITE ! 6 , 101 1 ) 

IF!  A! 25  > . GT .0 .  >  WR ITE 1 6 , 10 1 5 > 

PQRMAX- 1 . E+6 

IF!  A! 26 ) . GT . 1 .  )  PQRMAX -A! 26  >*A! 26 ) 

DUMMY-SORT! PQRMAX > 

PFLAG-. FALSE. 


QRFLAG-. FALSE. 

I F I NAL  =A ( 20 ) +  2 
WFINAL=A< IFINAL) 

WR I T  E  <  6 , 1 10)  A<33) .DUMMY. WF INAL 
110  FORMAT!-  RUN  WILL  STOP  AT  TIME  GT*.F5.1, 

1  '  TOTAL  BODY  RATE  GT’,F7.1.'  OR  WEIGHT  LT  -.F7.1) 
DELT-A! 30) 

DEL  PT=A ! 3 1 > 

TPR I  NT  =A ( 90 ) 

DDX-0. 

DDY=0. 

PT=0 . 

I  PTF=0 
I  PR  I  NT-99 
T IME  =0 . 

IT  I ME  =  -2 

DELT2=DELT*DELT/2. 

I E  =  —  1 
TG=A ! 33  > 

C 

COMMENT  -  CONTINUAL  CALCULATIONS 

2  CONTINUE 

IF  (TG.GT.DELT)  GO  TO  3 
DELT  =  TG 
I  E=  1 

3  ITIME-ITIME+1 

IF  ( I T I  ME  )  6,5,4 

4  I T I ME=  1 

5  TIME-TIME+DELT 

6  IF  ( T I  ME . GT . A ( 33 ) )  IE  =  2 
TSEQ=TIME-A(90) 

C 

CALL  KINE 
C 

CALL  AERO 
C 
C 

I F  <  T IME . LT . A ( 4  1  )  >  GO  TO  15 
C  BYPASS  VARIABLE  DELT 

TEMP=AP*AP  +AQ*AQ+AR  *AR 
IF(TEMP.LT. 10000. )TEMP=10000. 

DELT=10./SQRT(TEMP> 

DELT2=0.5*DELT*DELT 

15  IF!  <AP*AP+AQ*AQ+AR*AR ) .GT. PQRMAX  )  IE-3 

IF!  AW.LT.WFINAL  )  IE=4 
C 

GAMP-ATAN ! -AVZ/AVX ) *DPR 
GAMY=ASIN(AVY/AVT)*DPR 

TF  =  ( -AVZ ) /A ( 1 5 ) +SQRT ! < -AVZ > * < -AVZ > +  2  .  *A < I5)*AH)/A! 15) 

XF=AVX*TF+AX+DDX 

YF  =AVY*TF  +AY  +  DDY 

IF!  IPTF.NE.0  )  GO  TO  20 

I PTF  =  99 

XF0=XF  +0.001  +  A! 101 > 

YF0-YF+A! 102) 

DDX=A( 101 ) 

DDY=A! 102) 

20  DELTAX=!XF-XF0)*.3048 

DELTAY=!YF-YF0)*.3048 

D TOTAL -SORT !DELTAX*DELTAX+DELTAY*DELTAY> 

HAPO-AH+0. 5*AVZ*AVZ/A! 15) 

QRMX-SQRT I AQ*AQ+AR*AR ) 

IF  ( PFLAG )  GO  TO  52 
IF1AP.LT. A!28) >GO  TO  52 


PFLAG* . TRUE . 

WRITE (6, 51 )A(28) 

51  FORMAT 1 //5X , ’ ***ROLL  RATE  EXCEEDS F7 . 1  ) 
GO  TO  1002 

52  IF (QRFLAG )  GO  TO  54 
1F(QRMX.LT. A(29))GO  TO  54 
QRFLAG-. TRUE. 

WR I TE 1 6 , 53 )A( 29 ) 

53  FORMAT1//5X, ’***  QR  RATE  EXCEEDS ’,  F7  .  1  ) 
GO  TO  1002 

54  CONTINUE 

IF  ( At  3 1 ) .LE.0. )  GO  TO  9 
IF1TIME.GE.A1205) )  DELPT  =  At  32  > 

I F ( ABS( T1ME-TPR 1  NT ) . LT. DELT/2 .  )  GO  TO  22 

IF  (IE.GT.0)  PT-TIME 

IF  ( TIME . LT . PT  >  GO  TO  2 

PT-PT+DELPT 

IPOINT-1POINT+1 

1 F  t  1POINT.LT. 10)  GO  TO  1002 

IPO1NT-0 

WRITE(6. 1001 ) 


FORMAT!’  TIME’, 

,/9X, 

’psr 

’  ,9X,  ’ 

PHI 

’  ,9X. 

,  ’THETA’ , 

1 

9X ,  ’ 

X’ , 1 IX,  ’Y 

MIX 

.  ’Z’  , 

,  1  IX. ’ 

P’  . 

1 1 X  , 

1 Q  ’  ,  1 ).  X ,  ’  R  ’  , 

1 

4X ,  ’ 

Q  +  RIRSS)  ’  , 

2 

/10X 

, ’VX’ , 10X, 

!  ’ VY ’ . 

,  10X, 

,  ’VZ’ , 

10X 

.  ’AX’ 

1  , 10X, ’AY’ . 10X, ’AZ’ , 

3 

8X ,  ’ 

PDOT ’ ,8X, ’ 

1 QDOT 1 

1  ,8X, 

’ROOT 

’  ,/ 

10X,  ' 

■FI’ , 10X, ’ F2 ’ ,10X, 

4 

’  F3  ’ 

, 10X, ’Ml ’ , 

,  10X, ' 

’M2’  , 

,  10X .  ’ 

M3’ 

11X. 

’W’ ,9X, ’ 1XX’ ,9X, ’ 1 YY 

1 F  t  A1251.LT. 0.5  )  GO  TO  1002 
WR1TE16, 1011 ) 

101 1  FORMAT 19X, ’VXM’ ,  9X ,  ’ VYM' ,9X, ’VZM’  . 10X , ’XM ’ , 10X , ’ YM ’ , 

1  10X, ’ZM’ . 10X, ’TF ’ ,6X, ’DEL TAX' .  6X  , 'DELTAY' ) 

WRITE16, 1015) 

1015  FORMAT! 9X, ’ C 1 X ’ , 9X . ’ C 1 Y ’ , 9X . ’ C 1Z ’ . 8X , • D 1 VE ’ ,8X, ’TURN’ , 

1  8X , ’ HAPO ’ ) 

1002  WR1TE16.10I)  TIME 

101  FORMAT 1/2X,F10.6> 

WRITE  16. 102)  APS  1 , APH 1 . ATHETA , AX . AY , AZ . AP , AQ , AR , QRMX 
WR1TE16. 102)  AVX , AVY , AVZ , AAX , AAY , AAZ , APD , AQD , AR D 
WRITE  16. 102)  FI .F2.F3.M1 , M2 , M3 , AW , A1XX ,A1YY 

102  FORMAT! 10F12. 3) 

C  CUTOFF  TIME  FOR  TRANSFER  INTO  BALLISTIC  ROUTINE 

1 F 1  TIME  .GE.  A( 49 ) )  THEN 
VXS-5 . *  1 AC2X*AQ  ♦  AC3X*AR ) 

VYS-5 . *  1 AC2Y*AQ  ♦  AC3Y*AR  ) 

VZS-5 . *  1 AC2Z*A0  ♦  AC3Z*AR ) 

AVX -AVX  ♦  VXS+A142) 

AVY-AVY  ♦  VYS+Al 43 ) 

AVZ -AVZ  ♦  VZS+A144) 

A1 297 ) -8 . 8 

WR1TE16.331)  A142) ,A!43) ,A(44) 

331  FORMAT! IX, ’A142)  ’.F12.6,’  A 1 43 )  ’.F12.6,’  AA144)  ’, 

*  F12.6) 

WR1TE16.332)  VXS.VYS.VZS 

332  FORMATUX,’  VXS  ’.F12.6,’  VYS  ’.F12.6,’  VZS 

*  F12.6) 

WR1TE16.333)  AVX, AVY, AVZ 

333  FORMAT! IX, ’ DVX- ’ , F 1 2 . 6 , 5X , ’ DVY- ’ , F 1 2 . 6 , 5X , ’ DVZ= ’ ,F12.6 

CALL  BALLST1A) 

GOTO  515 
END  1 F 

C  CONVERT  TO  METRIC  FOR  OPTIONAL  PRINT 

IF!  A! 25 ) . LT . 0 . 5  )  GO  TO  22 
XM-AX* . 3048 
YM-AY* .3048 


'-'A’V 


jnyi  jp ;%  ,-y , 


2f “AZ* . 3048 
VXM*AVX* .3048 
VYM*AVY* .3048 
VZM*AVZ* .3048 
HAPO*HAPO* . 3048 

WRITE <6, 102)  VXM , VYM , VZM ,XM,YM,ZM,TF, DELTAX, DELTAY 
WR I TE ( 6 , 1020 )  AC  1 X , AC  I Y , AC  I Z , D I VE , TUR  N , HAP  0 
1020  FORMAT (3FI2.6.6FI2.3) 

C 

22  CONTINUE 
9  CONTINUE 

IF  ( IE.LT.0)  GO  TO  2 
ENDFILE  NT 
C 

I F {  A { 27 ) . GT . 0 . 5  )  WRITE<3,66)  DELTAX , DELTAY , DTOTAL 
66  FORMAT (/3X, 'DELTAX =’ ,FI2. 1  ,  '  DELTAY= ’ , F  1 2 . 2  , 

1  /3X, ’TOTAL  DELTA  IIP=’,FI2.I> 

C 

IF  (NMRC-1)  1,68.69 

68  CONTINUE 

69  CONTINUE 

I F (  A(40).LE.0.5  )  GO  TO  1 

I F  <  NMRC.LT.(NRC-LRC)  )  WR I TE <  4 , 4 444 >  DELTAX , DELTAY , ZERO 
I F  <  NMRC.GE. (NRC-LRC)  )  WR I TE (  4 , 4444  )  DEL  TAX , DELTAY , ZN I NE 
4444  FORMAT ( 3F 10.0  > 

C 

GO  TO  1 
END 

$PROG  AERO 

SUBROUTINE  AERO 

C  - 

REAL  IXX. IYY.M.Ml , M2 , M3 , M2A , M3A 
EXTERNAL  SIGN 

COMMON  /ALLT/  DELT , T IME , I TI ME , DELT2 

COMMON  /AKIN/  DR , ADP , AD Y , DF 1 , D F2 , DF 3 , DM  1 , DM2 , DM3 , WX , WY , WZ 
COMMON  /AKOT/AI , A2 , A3 , AX , AY . AZ , VX , VY , VZ , X , Y , Z , PD , QD , RD , P , Q , R , 

IPS  I .THETA, PH  I ,VA,VT,AMACH, ALPHA, BETA, ETA, QDP.H.W.CG, IXX, IYY.CA.CMT 

2, CLP,CLD,CLI,CMP,CMY,CNP,CNY,CIX,C2X,C3X,CIY,C2Y,C3Y,C1Z,C2Z,C3Z 

3, PHID,THED,PSID 
COMMON/DATA/  AI2800) 

COMMON  /CAERO/  CNA , CND , CMA , CMD , V23S , V23 , CAFAC , 

1  AVI , AV2 , AV3 , AV 4,ABSA,ABSB,SAD,SBD,XX,CL, 

2  DPS,DYS,DCY,DCG,D2V,V2A, V3A , SA , SB , F 1 , F 2 , F3 , M 1 ,M2,M3 
COMMON  /APLOT/  AL PHAA , BETAA , DPAERO , DYAERO , SADA , SBDA , 

&  CMPA.CMYA.CNPA.CNYA.PI , P 2 , P 3 , P 4 , P 5 , P6 , M2A , M3A 
COMMON  /CPLOT /  HERR , DFLARE , HGO , ADPX , AD YX , E LOS 
DATA  DPR/57.2957795/ 

C 

IF  ( ITIME.GE.0)  GO  TO  I 

C  CONVERT  FROM  INCHES  TO  FEET , +  OFFSET  GIVES  +  MOMENT 

DZCG=A ( I  3 ) / 1 2 . 

DYCG=A( 1 4 ) / 1 2 . 

W=A(  19) 

M=A< I9)/32. 17 

WX=0 

WY=0 

WZ=0 

ALPHA=0 . 

BETA=0. 

ETA=0. 

C  INPUT  ADDRESS  POINTERS  FOP.  FUNCTIONS  OF  WEIGHT  IN  LBS 

I XCG=A( 20  > 

I IXX=A<  2  I ) 

I IYY=A( 22  > 
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TT.  !» "T' 


ijl 

cl 


•V 


L 

\ 


*\ 


XGP  =A ( 23  ) 

SXOB=A( 38  ) 

TSTACK=A  <  34  ) 

C  INPUT  THE  IN-LB/PSI  COEFFS  FOR  ROLL  AND  PITCH/YAW 

RTR=A<  50) 

PYTR=A<  51 ) 

PHIPY=A(52)/DPR 

C  COMPUTE  COEFFS  TO  CONVEFiT  CHAMBER  PRESSURE  TO  TORQUE 

CEPT=COS(PHIPY)*PYTR/12. 

CEYT=SIN<PH1PY)*PYTR/12. 

CERT=RTR/12 . 

1  CONTINUE 

C 

COMMENT  -  6  DOF  ANGL ES -OF -ATTACK 
VXA=VX-WX 
VYA=VY-WY 
VZA=VZ-WZ 

VA=SQRT<VXA**2+VYA**2+VZA**2) 

V I A=C I X*  VXA+C 1Y*VYA+C 1Z  * VZA 
V2A=C2X*VXA+C2Y* VYA+C2Z*VZA 
V3A=C3X*VXA+C3Y*VYA+C3Z*VZA 
C 

c  *****  ROTATE  V2A , V3A , Q , R  INTO  AERO  AXES  •**•••«»* 

V2AS  =  V2A 
V3AS  =  V3A 

V2A  =  . 707*(  V2AS+V3AS  ) 

V3A  =  .  707*(  -V2AS  +  V3AS  ) 

OX  =  Q 
RX  -  R 

0  =  .  7J0T7 *  <  OX+RX ) 

R  *  .707M-OX  +  RX) 

I F  (  V2A*V3A  )  18J0T,  190,  180 
I  80  ALPHA=ATAN2<V3A, V1A)*DPR 
BETA =ATAN2<V2A, VIA) *DPR 
ETA=ATAN2  <  SORT! V2A**2+V3A**2 ) , VIA ) *DPR 
190  ALPHAA  =  ALPHA 
BETAA  =  BETA 
C 

I F (  TIME.GE.TSTACK  )  GO  TO  300 
FF=FN1V<A, 500, TIME) 

FR=FN1V<A, 550, TIME) 

WDOT«FNIV<A, 600, TIME) 

GO  TO  400 

300  TIME2»TIME-TSTACK 

FF«FN1V(A,I500,TIME2) 

FR-FNIVIA, 1550.T1ME2) 

FDR=FN1V(A, I700,TIME2  > 

WDOT-FN I V<  A , 1 600, TI ME2 ) 

C  UPDATE  MASS  PROPERTIES  HERE  INSTEAD  OF  IN  KINE 

400  I F  <  ITIME.LT. 0  )  GO  TO  555 

W=W-WDOT*DELT 
M*W/32 *  1 7 

C  INTERPOLATE  MASS  PROPERTIES 

555  XCG»FNIV<A,IXCG,W> 

IXX“FNIV<A, 1 1 XX , W > 

IYY»FN1V(A, I IYY.W) 

XBGP«<XGP-XCG)/I2. 

XB0B*(XCG-SXUB)/12. 

C  F20B-FNIVIA, 1650, TIME) 

C  F30B«FNIV<A, 1700, TIME > 

DP-FN1VIA, 305, TIME) 

DY-FNIVIA, 320, TIME  > 

DELTAT»SORT (DP*DP+DY*DY) 

F1N-C0S<DELTAT/DPR)*FF 


F2N=-SIN(DY/DPR)*FF 
F3N=-SIN(DP/DPR )*FF 
F1  =  T  1N-FR-FDR 
F2=F2N 
F3=F3N 

TRATI0=F1/A( 160) 

MI=A{ 1 6 1 ) *TRAT 10 
M2=A( 162) *TRAT 1 0 
M3=A ( 1 63 ) *TRAT10 
IF  <  TIME . LT.A( 164) )M2 =0. 

1 F ( TIME . GT . A( 165) >M2*0. 

C 

COMMENT  -  6  OOF  ACCELERATIONS 
A 1  =F 1 /M 
A2=F2/M 
A3=F3/M 

AX=C1X*A1+C2X*A2+C3X*A3 
AY=CIY*A1+C2Y*A2+C3Y*A3 
AZ=C1Z*A1+C2Z*A2+C3Z*A3+A( 15) 

Q  =  QX 
R  =  rx 

PD=DPR*M1 /  1XX 

QD  =  { M2*DPR  +  < 1YY-1XX)*P*R/DPR)/IYY 
RD= ( M3*DPR+ ( 1XX-1YY)*P*Q/DPR)/IYY 
V2A  =  V2AS 
V3A  =  V3AS 

1 F {  V2A*V3A  )  280.230,200 
280  ALPHA=ATAN2< V3A. VIA >*DPR 
BETA=ATAN2 ( V2A , VIA ) *DPR 
290  ADPX  =  . 707*<  DP -DY ) 

ADYX  =  . 707*( DP  +  DY  > 

RETURN 

END 

SPROG  SOLVE 

SU8P0UT 1 NE  SOLVE(X,Y,XIN,XOUT , A , B , C ) 

DIMENSION  X <  3 ) .Y<3> 

C  COMPUTE  QUADRATIC  COEFFS  FOR  Y=A*X**2  +  B*X  +  C 

D1=Y< 1 )/( (X( 1  )-X(2)  )*(X( 1 >-X(3>  )  ) 

D2  =  Y(2)/( (X(2)-X( 1  )  )*(X<  2  )-X<3>  >  > 

D3  =  Y(3)/{ (X!3)-X< 1  )  )  * < X ( 3  ) -X < 2 >  >  ) 

C  =  D1*X(2)*X(3)+D2*X( 1 > *X < 3 ) +D3*X (  1  )*X(2> 

B  =  -D1*<X<2)+X<  3) )-D2*<X( 1 )+X<3 ) > -D3* < X {  1  )  +  X( 2)  ) 
A=D 1 +D2+D3 

X0UT=A*XIN*X1N+B*X1N+C 

RETURN 

END 

SPROG  FN1V 

FUNCTION  FNIV  <A,I,X> 

C  - 

COMMENT  -  FUNCTION  OF  ONE  VARIABLE  -  4/15/70 
DIMENSION  A ( 2800  > 

N*A { 1 )+.5 

IF  (N.EQ.0)  GO  TO  4 
IF  (N.EO. 1 >  GO  TO  5 
DO  1  0=4, N, 2 
K=  I  +0 

F N 1 V=A <  K+ I  ) 

DEL=X-A(K) 

IF  (DEL)  2,3,1 

1  CONTINUE 

2  CONTINUE 

IF (A( K) . EQ.A( K-2  > )  GO  TO  3 

SLOPE” (  FNIV  —A ( K— I ) )/(A(K)-A(K-2) > 

FNIV=FNIV+DEL*SLOPE 
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3  return 

4  FN1V-0. 

RETURN 

5  FNlV-Ad  +  2) 

RETURN 

END 

SPROG  INPUT 

SUBROUTINE  INPUT  (A, NDIM, IN1TFN , NR > 

C  - - - 

COMMENT  -  DATA  INPUT  SUBROUTINE  -  1/2B/7B 

DIMENSION  LMC  (  200)  ,  NO  IST(  200)  ,  VAL  1  { 200)  ,  VAL2  (  200) ,  VAL0<  200)  ,  IRC<  2  ) 
DIMENSION  A(NDIM) , ISEQ< 100) ,RSEO( 100 ) , IS < 3 ) , NS ( 3 ) , VSE ( 5 , 6 ) 
EQUIVALENCE  ( ISEQ< 1 ) ,RSEQ< 1 ) ) 

LOGICAL  PSEQ , POI NT , AZERO , DA , LA , CARLO , PMC 
COMMON  /MONT/NMRC , LRC ,NRC , IB 
I NTEGER*2  I N < 72 ) , SYMBOL < 45 ) , I T , LABEL ( 72 > 

DATA  SYMBOL  /1H  , 1 H0, 1 H 1 , 1 H2 . 1 H3 , 1 H4 , I H5 , 1 H6 , I H7 , I HB , I H9 , I H . , 1 H+ , 
11H-,IH*,1H/,1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,IHI,1H0,1HK,1HL,IHM, 
21HN.1H0, 1HP, 1HQ, 1HR, 1HS, 1HT, 1HU, 1HV, 1HW, 1HX, 1HY, 1HZ, 1H( , 1H> , 

3  1H’/ 

DATA  0SEQ,NSEQ,LSEQ,PSEQ,A2ER0/£r,£r,  1 ,2*  .FALSE./ 

DATA  DA, CARLO, PMC, MRC.NMC, II ,MMM/3*. FALSE. ,4*0/ 

DATA  VSE/- 1 . ,-.6, -.2,. 2,. 6,. 2, -.6,-1.,. 6, -.2,. 6,-1,, -.2, -.6,. 2, 

1  .6,— .2,— .6, .2,— 1.,— .6,— .2, .2,— 1.,. 6,— 1.,. 2, .6,— .2,— .6/ 

C 

COMMENT  -  INITIALIZE 

IF  (AZERO)  GO  TO  2 

INFILE-5 

NMRC-0 

LRC-J0T 

NRC-fiT 

NR-0 

DO  1  1-1, NDIM 

1  A< I )«0. 

AZERO-. TRUE. 

C 

COMMENT  -  IDENTIFY  AND  SET  UP  COMMANDS 

2  I F ( CARLO )  GO  TO  500 

IF  (NSEQ.GT.jBT)  GO  TO  450 
WRITE (6, 4) 

4  FORMAT  (5HJ0TDATA) 

3  R EAD (INFILE, 10)  IN 

10  FORMAT  ( 72A1 ) 

WRITE (6, 11)  IN 

11  FORMAT  (1X.72A1) 

0-1 

NB-0 

IB  IT-IN(O) 

DO  20  1-3,11 

IF  ( I T. EQ. SYMBOL ( I ) )  GO  TO  100 
20  CONTINUE 

C  ALPHABETIC  CONTROL  STATEMENTS 

IF  ( IT. EQ. SYMBOL (34))  GO  TO  200 
IF  ( IT.EQ.SYMBOL(2) )  GO  TO  50 
IF  ( IT.EQ.SYMBOL(35) )  GO  TO  70 
IF  (IT, EQ. SYMBOL ( 1 ) )  GO  TO  B0 
IF  ( IT.EQ.SYMBOL(29))  GO  TO  90 
IF  ( IT. EQ. SYMBOL! 15) )  GO  TO  40 
IF  ( IT. EQ. SYMBOL! 19))  THEN 
CALL  BALLST(A) 

IB-1 
RETURN 
END  IF 

IF  ( IT.EQ.SYMBOL(21 ))  STOP 


102  CONTINUE 

IF  (LOC.GT.NDIM)  GO  TO  30 
IF  < I N  <  J  > . EQ.SYM80L! I ) )  GO  TO  115 
IF  < IN  <  J ) . EQ .SYMBOL  <  29  > )  GO  TO  109 
IF  < IN!0).EQ.SYMBOL<35> )  GO  TO  112 
C 

COMMENT  -  SET  UP  FUNCTION  MODIFICATION 
IF  <  LOC . LT . IN ITFN )  GO  TO  30 
N  1  -A ( LOC ) 

N2«A<  LOC+I ) 

NT-N1+N2 

NX-NI+N2 

IF  <  NT . LE .0)  GO  TO  30 

IF  !IN!0).NE.SYMBOL!22))  GO  TO  107 

L-3 

IF  <  N 1  . EQ . 1 >  L-2 
IF  (NX.GT.0)  L-NI+3 
L  I  "2 

IF  (NX.GT.0)  L I =-N 1 - 1 
LASTL-NX+NT+1 
GO  TO  120 

107  IF  < I N<  J  ) . NE . SYMBOL  <  25 ) )  GO  TO  108 
IF  <  N 1 . LT . 2 )  GO  TO  30 

L-2 

LI-2 

IF  (N2.GT.0)  LI-1 
LASTL-N1+2-LI 
GO  TO  120 

108  IF  < I N ( 0 ) .NE. SYMBOL! 18) )  GO  TO  30 
IF  (N2.LT.2)  GO  TO  30 

L-NI+2 

LI-NI+I 

IF  (N1.EQ.0)  LI-2 
LASTL-NX+NT+2-L I 
GO  TO  120 
C 

COMMENT  -  SET  UP  MONTE  CARLO  DATA 

109  IF  (LOC.GT. INITFN)  GO  TO  30 
O-O  +  I 

DO  110  1-3,5 

IF  (IN(O).EQ. SYMBOL  < I  )  )  GO  TO  III 

110  CONTINUE 
GO  TO  30 

111  0-0+1 

IF  < IN<0 ) .NE. SYMBOL!  I )  )  GO  TO  30 
NMC-NMC+1 

IF  (NMC.GT.200)  GO  TO  30 
LMC ! NMC ) -LOC 
NDIST!NMC)»I-2 
MMM-I 
GO  TO  1 15 
C 

COMMENT  -  SET  UP  SEQUENCED  DATA 

112  ISEQ(LSEQ)-LOC 
NSEQ-rSEQ+1 
OSEO-LSEQ+I 
LSEQ-LSFQ+2 

C 

COMMENT  -  SET  UP  TO  READ  DATA 
115  L-0 
0P--1 
LASTL--1 
120  NB-I 


•l 


LASTOP-1 
122  NO-0 
NP=0 
FN=0. 

POINT-. FALSE , 
NI-0 


COMMENT  -  READ  DATA  AND  COMPUTE  VALUES 
125  0*0+1 

IF  ( IN ( 0 ) . EQ .SYMBOL  < I ) )  GO  TO  150 
NB-0 

DO  135  1=2,11 

IF  (IN(J).EQ. SYMBOL ( I ) )  GO  TO  140 
135  CONTINUE 
GO  TO  160 
140  N I = 1 

NO=I0*NO+ 1-2 
IF  (POINT)  NP-NP+I 
145  IF  (0.LT.72)  GO  TO  125 
150  NB-NB+1 

IF  (NI)  182,152,180 
152  IF ( NB . LT. 12)  GO  TO  145 

IF  (L.GT.LASTL)  GO  TO  187 
155  READ  ( INFILE, 10)  IN 
0=0 
OP--1 
GO  TO  120 

160  IF  (IN(J).NE. SYMBOL ( 1 2 ) >  GO  TO  165 
POINT-. TRUE. 

GO  TO  145 

165  IF  (IN(O).NE. SYMBOL ( 38  ) )  GO  TO  170 
POINT*. FALSE. 

LV-LOC+L 

FN-A(LV) 

NI--1 
GO  TO  145 

170  DO  171  1-13,16 

IF  (IN(O).EQ. SYMBOL ( 1 ) )  GO  TO  179 

171  CONTINUE 

IF  <IN(J).NE.SYMBOL<43)>  GO  TO  172 

OR-O 

OP -NO- I 

IF  (LOC.GE. INITFN.AND.NO. EQ.0)  OP-9999 
IF  (0.EQ.72)  GO  TO  30 
GO  TO  120 

172  IF  (IN(O).NE .SYMBOL ( 44 ) >  GO  TO  30 
IF  (OP)  30,150,174 

174  OP-OP-1 
O-OR 

GO  TO  150 

179  NEXTOP -I - 1 2 

IF  (NI.LT.0)  GO  TO  182 

180  FN-NO 

F  N-FN/ 10. **NP 

182  IF  (LASTOP .EQ. I )  VAL-VAL+FN 
IF  (LASTOP. EQ. 2)  VAL-VAL-FN 
IF  (LASTOP. EQ. 3)  VAL-VAL  *FN 
IF  (LASTOP. EQ. 4)  VAL-VAL/FN 
IF  (NB.EQ. 1 )  GO  TO  183 
LASTOP-NEXTOP 
IF  (0.EQ.72)  GO  TO  30 
GO  TO  122 
C 

COMMENT  -  STORE  SEQUENCED  VALUES 


183  IF  (OSEO.EQ.0)  GO  TO  188 
IF  (LSEQ.GT.100)  GO  TO  186 
RSEQ<LSEQ)-VAL 

IF  (LSEQ.EQ. JSEQ+1 )  1SEQ( LSEQJ-VAL 

LSEQ-LSEQ+1 

GO  TO  120 

186  OSEQ-0 
LSEQ-OSEQ-1 
NSEQ-NSEQ-1 
GO  TO  30 

187  IF  <3SEQ.EQ.0>  GO  TO  3 

IF  ( LSEQ . LT . OSEQ+4 )  GO  TO  186 
lSEQ{0SEQ)-LSEQ-0SEQ-2 
OSEO-0 
GO  TO  3 
C 

COMMENT  -  STORE  MONTE  CARLO  VALUES 

188  IF  (MMM.NE. 1 )  GO  TO  192 
11-11+1 

IF  (ll.GT.l)  GO  TO  190 
VAL1 (NMC)-VAL 
GO  TO  120 
190  VAL2(NMC)-VAL 
1  1-0 
MMM-0 
GO  TO  3 
C 

COMMENT  -  STORE  DATA  VALUES 

192  LV-LOC+L 

IF  (LV.GT.ND1M)  GO  TO  30 
A<  LV ) -VAL 

IF  (LOC.LT. 1N1TFN)  GO  TO  193 
IF  (L-l)  196,194,197 

193  IF  (0-72)  196,3,30 

194  Ll-1 

LASTL-(A( LOC )+ 1 . )*(A(LOC+l  )*1 .  ) 

196  L-L  +  l 

GO  TO  120 

197  IF  (L.GE.LASTL)  GO  TO  3 
IF  (L1.LT.0)  GO  TO  198 
L-L+Ll 

GO  TO  199 

198  IF  ( MODI LASTL-L ,-Ll>.EQ.0)  L-L+l 
L-L+l 

199  IF  (0-72)  120,155,30 
C 

COMMENT  -  SET  RUN  NUMBER.  PRINT  VALUES,  AND  START  RUN 

200  NP-0 
NB-0 

LA-. FALSE. 

210  0-0+1 

IF  ( 1N(0) . EQ. SYMBOL ( 1 ) )  GO  TO  270 
NB-0 

DO  230  1-2,11 

IF  ( 1 N( 0 ) . EQ. SYMBOL ( 1 ) )  GO  TO  250 
230  CONTINUE 

IF  ( IN(0).NE.SYMBOL(28>)  GO  TO  240 
LA- .TRUE . 

NA-0 

GO  TO  210 

240  IF  ( 1N(O).NE.SYMBOL(20))  GO  TO  30 
DA-. FALSE. 

GO  TO  210 

250  IF ( LA )  GO  TO  265 
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NP=I0*NP+I-2 
GO  TO  210 
265  NA»I0*NA+I-2 
GO  TO  210 
270  NB=NB+ 1 

IF  (NB.LT.I2)  GO  TO  210 
I F  <  LA )  DA  =  LA 
I F  <  NA . EQ . 0)  NA-NOIM 
IF  (NP.GT.0)  GO  TO  278 
NR=NR+I 
GO  TO  280 
278  NR-NP 

280  IF  ( NSEQ.GT.0)  GO  TO  400 
282  IF  (CARLO)  GO  TO  500 

291  IF< .NOT. DA)  GO  TO  295 
0=0 

K=9 

IF  (K.GT.NA)  K  =  NA 

WRITE  (6,292)  0 , NR , ( A( I ) , I • I , K ) 

292  FORMAT  (I2H  INPUT  ARRAY//I6.5H  RUN , 1 6 , 2X , 9F 1 1 . 4 ) 

293  IF  (K.EQ.NA)  GO  TO  295 
0=0+10 

K=K+I0 

IF  (K.GT.NA)  K=NA 

WRITE  (6,294)  0 , ( A( I ) , I =0 , K > 

294  FORMAT  ( 1 6 , 2X , I0F 1 1 . 4 ) 

GO  TO  293 

295  IF  (NSEQ.EQ.0)  RETURN 
IF  ( . NOT . PSEQ  )  RETURN 
1  =  1 

0=1 

WRITE  (6,297) 

297  FORMAT  (I7H  SEQUENCEO  VALUES) 

298  L  =  I SEQ ( I ) 

WRITE  (6,299)  L.A(L) 

299  FORMAT  ( 1 6 , 2X , I0F I  I . 4/ ( 6X , 10F 1 1 . 4 ) ) 

I = I + ISEQ( I+I )+3 

0=0  +  1 

IF  (O.GT.NSEQ)  RETURN 
GO  TO  298 
C 

COMMENT  -  INITIALIZE  SEQUENCED  VALUES 
400  DO  405  1=1,3 
IS (  I  )  =  I 
405  NS( I  )“0 
1  =  1 
0=1 

410  L-ISEQ(I) 

A(L )=RSEQ( 1+3) 

K=ISEQ( 1+2) 

NS ( K )  =  ISEQ ( I  +  I > - I 
0=0  +  1 

I  =  I +NS ( K ) +4 

IF  (O.GT.NSEQ)  GO  TO  282 
GO  TO  410 
C 

COMMENT  -  UPDATE  SEQUENCED  VALUES 
450  DO  455  1=1,3 
IS( I ) = IS ( I )+I 

IF  ( IS ( I ) . LE  .  NS( I  )  )  GO  TO  465 
455  IS ( I ) = I 
I»I 
0=1 

IF  (PSEQ)  WRITE  (6,297) 
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L=ISEQ< 1 ) 

1  =  1  +  ISEQ( 1+11+3 
A(L  >=RSEQ( 1-1 > 

0=0  +  1 

IF  (PSEQ)  WRITE  (6,299)  L.A(L) 
IF  (O.LE.NSEQ)  GO  TO  460 
NSEQ=0 
LSEQ=1 
OSEQ=0 
GO  TO  3 
CONTINUE 
NR=NR+ 1 
1  =  1 
0  =  1 

K= 1SEQ ( I +2 ) 

K=lS(K)+I+2 
L  =  1  SEQ( 1 > 

A(L  )=RSEQ(K) 

1=I+ISEQ( 1+1 >+3 
0=0  +  1 

IF  (0.GT.NSEQ)  GO  TO  282 
GO  TO  470 


VALUES 
TO  530 


504 


COMMENT  -  UPDATE  MONTE  CARLO 
500  IF  ( LRC+MRC . GE . NRC )  GO 
MRC=MRC+1 
NMRC=MRC 

CALL  RANI (MRC+LRC ) 

NC=0 

IF  (.NOT. PMC)  GO  TO 
WR I TE ( 6 , 502 ) 

502  FORMAT  (//18H  ARRAY  LOCAT I  ON , 5X , 1 3HNOMI NAL  VALUE, 
1  5X.9HRUN  VALUE ,5X, 10H DIFFERENCE) 

504  NC-NC+1 

VI =VAL 1 ( NC ) 

V2=VAL2 ( NC ) 

ND=ND 1ST ( NC ) 

LL-LMC(NC) 

IF  (MRC.EQ.l)  VAL0( NC )»A( LL ) 

NZ=MRC+LRC 

NSE=MOD( (NZ-1 )/5+NC-l , 120) 

1SE=M0D (NSE,6)+1 
OSE=MOD ( NSE , 5 ) 

KSE=MOD ( NSE+NSE/60 , 4 ) + 1 
LSE=MOD(OSE+(NZ-1 )"KSE,5)+1 
RUN  IF  =RANU ( . 4 ) +VSE ( LSE , 1SE ) 

IF  (ND.NE. 1 )  GO  TO  506 

A(LL )=V1+V2*RUNIF*( 1 . 1075+ . 1 42/ ( 1 .0898-ABS ( RUNI F ) ) ) 
GO  TO  520 

IF  (ND.NE. 2)  GO  TO  510 
A(LL )=V1+V2*RUN1F 
GO  TO  520 

IF  (RUNIF.LT. 0. )  GO  TO  512 
A(LL )=V2 
GO  TO  520 
A(LL)»V1 

IF  ( .NOT. PMC)  GO  TO 
DIFF=A(LL  >-VAL0(NC) 

WRITE (6, 522)  LL , VAL0(NC ) ,A(LL ) , DIFF 
FORMAT  (8X.14, 1 2X , F 10 . 3 , 6X , F 10 . 3 , 4X , F 10 . 3 ) 

IF  (NC.LT.NMC)  GO  TO  504 
IF  (MRC.EQ.l  )  GO  TO  291 
RETURN 


506 


510 


512 

520 


526 


522 

526 
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COMMENT  -  RESTORE  MONTE  CARLO  VALUES 

530  DO  531  NC-l.NMC 
L  L  =  LMC  <  NC  > 

531  A<LL)-VAL0(NC> 

NMRC-0 

MRC=0 

IF  (NSEO.GT.0)  GO  TO  450 
CARLO-. FALSE. 

NMC-0 
GO  TO  3 
END 

SPROG  KINE 

SUBROUTINE  KINE 
C  - 

real  IXX  I yy 

COMMON  /ALLT/  DELT , TI ME . ITIME , DELT2 

COMMON  /AKIN/  DR . DP , DY . DF 1 . DF2 , DF3 , DMI , DM2 , DM3 , WX , VY, VZ 
COMMON  /AKOT/AI , A2 , A3 , AX . AY , AZ . VX , VY , VZ , X , Y , Z , PD , QD , RD , P , Q , R , 

I  PS  I .THETA, PHI , VA,VT,AMACH, AL PHA, BETA, ETA, QDP,H,V,CG, IXX , IYY.CA.CMT 

2,  CLP .CLD.CLI ,CMP , CMY , CNP , CN Y , C 1 X , C2X , C 3X , C 1 Y , C2Y, C3Y, C 1Z , C2Z , C3Z 

3, PHID,THED,PSID 
COMMON/DATA/  A(2800) 

DIMENSION  QA ( 4 ) , QAP ( 4 ) 

DATA  DPR/57.2957795/ 

IF  (ITIME. GE.0)  GO  TO  I 
X-A(l) 

Y-A( 2 ) 

Z--A(3) 

V-A( 4 ) 

CA5-COS( A( 5 ) /DPR ) 

VX  =V*CA5*COS(A<  G )/DPR ) 

VY  =V*SIN(A(5)/DPR) 

VZ  — V«SIN(A(6)/DPR)*CA5 
PSI -A( 7 ) 

THETA-A!  8 ) 

PHI-A( 9 ) 

CPSI-COS(A(7)/DPR) 

SPSI-SIN(A(7)/DPR) 

CPHI-COS(A(9)/DPR) 

SPHI-SIN(A(9)/DPR) 

CTHETA-COS(THETA/DPR) 

STHETA-S I N ( THETA/DPR ) 

CIX-CTHETA*CPSI 

CIY-SPSI 

C i2--STHETA*CPS I 

C2X— CPHI*SPSI*CTHETA+SPHI*STHETA 
C2Y-CPHI*CPSI 

C2Z -S PHI *C THETA ♦STHETA*SPSI*CPHI 
C3X-SPHI *SPSI*CTHETA*CPHI*STHETA 
C3Y— SPHI*CPSI 

C3Z“CPHI*CTHETA-STHETA*SPSI *SPH I 
IQUAT-0 

QAP ( 4 ) -0. 5*SQRT (I.+C1X+C2Y+C3Z) 

QAP4S-0. 25* ( I . +CIX+C2Y+C3Z ) 

880  FORMAT! IX, ’QAP 3’ > 

QAP(3)-SQRT(QAP4S-0.5*(CIX«-C2Y)  ) 

IF ( C2X .GT. C IY  >QAP ( 3 I--QAP ( 3  > 

QAP ( 2 )-SQRT( QAP  4S-0. 5*(CIX+C3Z) ) 

IF( CIZ.GT.C3X  >QAP ( 2 ) --QAP ( 2 ) 

QAP ( I )«SQRT(QAP4S-0.5*(C2Y+C3Z> ) 

IF (C3Y.GT.C2Z )QAP( I I--QAP ( I ) 

C 

P-A( 10) 

Q-A( II  ) 


R=A( 12) 

GO  TO  2 
C 

1  CONTINUE 

VX  -  VX  +AX*DELT 
VY  -  VY  +AY*DELT 
V2  «  V2  *A2*DELT 
X=X*VX*DELT*AX*DELT2 
Y=Y*VY*DELT*AY*DELT2 
2=2+V2*DELT+A2*DELT2 

p=p*pd*delt 

Q=Q*QD*DELT 

R=RfRD*DELT 

C 

DELP-(P*DELT+PD*DELT2)/DPR 
DELQ«(Q*DELT*QO*DELT2)/DPR 
DELR«(R*DELT+RD*DELT2)/DPR 
DELAS»DELP*DELP*OELQ*DELQ+OELR*DELR 
DQA-0. 1 25*DELAS* ( DELAS/48 . -  1 . )*1 . 

SQA-0.  5-DELAS/48  . 

HXA-DELP*SQA 

HYA-DELQ*SQA 

H2A-DELR*SQA 

QA( 1  )  =DQA*QAP ( 1 ) +H2A*QAP ( 2  )-HYA*QAP ( 3 ) +HXA*QAP ( 4 ) 
QA( 2 ) --H2A-QAP ( 1 ) +DQA*QAP ( 2 ) +HXA*QAP ( 3 >  +HYA*QAP ( 4 ) 
QA( 3 ) -HYA*QAP ( 1 ) -HXA*QAP ( 2 )+DQA*QAP ( 3  >*H2A*QAP ( 4 ) 
QA( 4 ) «-HXA*QAP ( 1 ) -HYA*QAP  <  2 ) -H2A*QAP ( 3  )+DQA*QAP ( 4 ) 
DO  444  1-1,4 

444  QAP ( I )»QA(  I  ) 

IQUAT - I QUAT ♦ 1 

I F  {  IQUAT.LT.  lJ0fj0f)GO  TO  446 

TEMP«j0T.5*{3.-QA(  1  )*QA(  1  ) -QA<  2  )  *QA  (  2  > -QA  (  3  >  *QA  (  3  )  - 
*  QA( 4  )*QA( 4 ) ) 

DO  445  1-1,4 

445  QAP ( I ) -TEMP*QAP ( I ) 

IOUAT-0 

446  CONTINUE 

C 1 X-QAP ( 1 ) **2-QAP (2)**2-QAP(3 ) **2+QAP ( 4 ) **2 
C 1 Y-2 . * ( QAP ( 1 ) *QAP  <  2 ) +QAP ( 3 ) *OAP ( 4 ) ) 

C 12-2 . * ( QAP ( 1 ) *QAP ( 3 ) -QAP {  2 ) *C AP {  4) ) 

C2X-2  .  *  ( QAP  (  1  )*QAP  (  2 )  -QAP  <  3  )  *Q>'P  (  4  )  ) 

C2Y-QAP {2)**2  +  QAP(4)**2-QAP{  1  )••<.  -QAP(3>**2 
C22-2 . * ( QAP { 2 ) *QAP ( 3 ) *QAP ( 1 ) *QAP ( 4 ) ) 

C3X-2 . * ( QAP ( 1 ) *QAP ( 3  )  +QAP ( 2 ) *QAP  <  4) ) 

C3Y=2.*(QAP(2) *QAP ( 3 ) -QAP ( 1 )*QAP<4  >  ) 

C32-QAP ( 3 ) **2*QAP (4)**2-QAP( 1 )**2-QAP(2)**2 
THETA-ATAN2(-C12,C1X)*DPR 
PHI-ATAN2(-C3Y,C2Y)*DPR 
PSI-ASIN(C1Y)*DPR 
C 

2  CONTINUE 
H--2 

VT  -  SORT ( VX**2*VY**2+V2**2 ) 

RETURN 

END 

SPROG  SIGN 

FUNCTION  SIGNIAl ,A2> 

Q-ABS ( A1 ) 

iF(Q.NE.jer.jer)GO  to  i 
SIGN-A1 
RETURN 
1  CONTINUE 
2-1  .  J0T 

IF(A2.LT.ff.0>2--l  .0 
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SIGN-Z*ABS<AI  ) 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SPROG  RANI 

SUBROUTINE  RANK  I  ,  N  ) 

C  INCLUDED  IN  ORIGINAL  PROGRAM  BUT  NOT  ACTIVATED  BY  UAH 

C  RANDOM  NUMBER  INITIATOR  FOR  8/32 

C  INITIALIZES  A  DOUBLE  PRECISION  FRACTION ,N ,  FOR  SUBSEQUENT  USE 
C  BY  THE  RANU  AND/OR  RANG  FUNCTIONS 
C 

DOUBLE  PRECISION  M, N 
M-I7I79B69184.D0 
C  MODULUS  M  =  2**34 

N-  <  56295  1 4 1 3 . D0  +  2B 1 82B 1 72 . D0* I ) /M 
N  =  N- I D I  NT  <  N ) 

C 

C  N  IS  NOW  A  DOUBLE  PRECISION  FRACTION 

C  N-2**<-34)  +  MISC  HIGHER  BITS  WHICH  DEPEND  ON  I 

C 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SPROG  RANU 

FUNCTION  RANU  <  DP FRC , A ) 

DOUBLE  PRECISION  DPFRC.XX 

XX-DPFRCM3I075.D0 

N=XX-IDINT(XX) 

Xl-N 

RANU-A*XI 

RETURN 

END 

ccccccccccccccccccccccccccccccccccccccccccccccccccc 

SPROG  BALLST 

SUBROUTINE  BALLST(A) 

DIMENSION  A<2800) 

COMMON  /ALLT/  DELT , TI ME , ITI ME , D ELT2 

COMMON  /AKOT/AI , A2 . A3 , AX . AY , AZ . VX , VY , VZ , X , Y , Z , P D ,QD , RD . P , Q, R , 

I  PS  I , THETA , PHI , VA , VT, AMACH , AL PHA . BETA , ETA . QDP , H , W, CG , I XX . I YY , CA, CMT 

2, CLP,CLD,CLI ,CMP, CMY ,CNP,CNY,CIX,C2X,C3X,CIY,C2Y,C3Y,C1Z,C2Z,C3Z 

3, PHID,THED,PSID 
BC-A<45> 

BDT-AI 46 ) 

BDP-AI 47  ) 

IPT-0 

AXO-AX  *DELT/BDT 

AYO-AY  *DELT/BDT 

AZO-AZ  *DELT/BDT 

VXO=VX 

VYO-VY 

VZO-VZ 

TNXTPR=TIME+BDT 

BDTI-BDT/2. 

WR ITE  <  6 , 707 )T IME ,X , Y ,Z , VX , VY , VZ , AX , AY.AZ , RHO 
1  VT=SQRT< VX**2  +  VY**2  ♦  VZ**2) 

CALL  ATMO ( -Z , RHO , VS ) 

AX--16. I*RHO*VT*VX/BC 
AY--I6. I*RHO*VT*VY/BC 
AZ— 16.  I*RHO*VT*VZ/BC  +  32.2 
VX-  VX+ (AX+AXO ) *BDTI 
IF<ABS(VX).LE.l.E-35)  VX-0.0 
AXO-AX 

VY-  VY+  < AY+AYO ) *BDTI 
IF(ABS<VY) .LE. I .E-35)  VY-0.0 


AYO=AY 

VZ  =  VZ+! AZ+AZO ) *BDT I 
AZO=AZ 

X  =  X+I VX+VXO)*BDTI 
VXO=VX 

Y=  Y  + ! VY+VYO ) *BDT I 
VYO=VY 

Z=  Z+!VZ+VZO)*BDTI 
VZO=VZ 

TIME=TIME+BDT 
IFUPT.EQ.  1  )  GOTO  767 
WR  I TE  <  6 , 727  ) 

727  FORMAT  (6X,  ’  T I  ME  ’  ,  9X  ,  ’  X  ’  ,  1 2X  ,  '  Y*  ,  9X  ,  ’  Z  \  1 2X  ,  ’  VX  ’  .  1J0TX  ,  ’  VY  ’  , 

+  1 IX, ’VZ* , 1 IX, ’AX* ,9X, ’AY'  ,9X, ’AZ’ ,9X, 'RHO*  > 

I P  T=  I 

767  IFCTIME  .GE.  TNXTPR )  THEN 

WR ITE ( 6 , 707 )T IME ,X , Y, Z , VX , VY, VZ , AX , AY, AZ , RHO 
707  F  OR  MAT  !10F12.3,F12.9) 

TNXTPR=TNXTPR+BDP 

ENDIF 

IF  (Z  .GE.  0.0)  THEN 

WRITE(6,757)TIME,X,Y,Z,VX , VY , VZ , AX , AY . AZ 
757  FORMAT! 10F12. 3) 

RETURN 
ENDIF 
GO  TO  I 
RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SPROG  ATMO 

SUBROUTINE  ATMO! H , RHO , VS  ) 

C  H  IS  MEASURED  POSITIVE  UP 

C  H  HAS  UNITS  OF  FEET 

C  REFERENCE  AEROPLANE  AERODYNAM ICS , DOMMASC H , 1 95 1 

C  HNA=HE I GHT  ABOVE  WHICH  NO  SENSIBLE  ATMOSPHERE  EXIST 

HNA  =  3j0Tj0j0j0j0.j0 
RHOSL-.j0Tj0T2378 

C  FOR  HNA  =  3000J0  FEET  THE  CENSITY  IS  ONE  MILLIONTH  THAT  OF  SL 

C  RHOSL  IS  ATMOSPHERIC  DENSITY  OF  SEA  LEVEL  IN  SLUGS/FT**3 

IFIH.LE. 35332.0)  THEN 

RR=! 1 .0-.0000068*H) **4.256 
ELSE  IFIH.LE.HNA)  THEN 

RR  =  I .32/EXP!  1 .452  +  ! IH-35332. J/20950. ) ) 

ELSE 

RR=0.0 
ENDIF.  ... 

t)Hrt-Rfi*RWnsl 

C  VS  IS  SPEED  OF  SOUND  IN  FT/SEC 

C  I F  t  Hv  t  E . ^&^^2  40 )  THEN 

C  VSTin6.4  -  . 0042*H 

C  ELSE  IFIH.LE. 83000.0)  THEN 

C  VS 4 9 B 8. 5 

C  ELSE  I F I H I LE . 173885.0)  THEN 

C  VS-968.5  +  . 00 1 5 1 5  * ! H - 8  3000 . ) 

C  ELSE  IF(H.LE.HNA)  THEN 

C  VS" 1 106 . 19-  .0025324*! H-I73885. ) 

C  ELSE 

C  VS-786.8 

C  ENDIF 

RETURN 

M.M.  HALLUM  13  JAN  1983  -  AERO  EQUATIONS 
END 


C 


CORRECT  EQUATIONS 
FOR  AERODYNAMIC 
AND  ATMOSPHERIC 
CONDITIONS 


TRANSFORM  MISSILE 
VELOCITIES  TO 
GROUND  FRAME 


CALCULATE 

DERIVATIVES 


TI^.GT.TPR 


DETERMINE  AERODY¬ 
NAMIC  AND  ATMOS¬ 
PHERIC  CONDITIONS 


^^IF^v 
^  (TIME.GT 
JTFINAD.OR. 
\jMPACT  ^ 


INTEGRATE 

VARIABLE 

DERIVATIVES 


PRINT  INTERMEDI¬ 
ATE  RESULTS 


PRINT  FINAL 
VARIABLE  VALUES 


2.7S  DIGITAL  SIMULATION 


S9ATCH 

SXREF 

SPROG  MAIN275 
C 


WRHD1S1 


PROGRAM  MAIN<INPUT/0UTPUT/TAPE5=INPUT/TAPE6=0UTPUT/TAPE2/TAPE10> 
***************************************** 

*  6-DOF  SIMULATION  * 

*  2.75  INCH  ROCKET  * 

*  AND  * 

*  DOWNWASH  PROFILE  * 

*  * 
*************  **************************** 

LOGICAL  VARIABLE  IMPACT  STOPS  TRAJECTORY  WHEN  ROUND  REACHES 
6R0UND/THAT  IS  ALTITUDE  XX<2/9>  GOES  POSTIVE. 

LOGICAL  VARIABLE  *  NOTIME5  STOPS  TRAJECTORY  WHEN  FLIGHT  TIME  HAS 
EXCEEDED  ALLOTED  RUNTIME  TIME  OF  *  TF INAL ' . 


LOGICAL  BURNOUT 

LOGICAL  LNEXIT 

LOGICAL  SET90 

REAL  IXO 

REAL  IX 

REAL  IXDOT 

REAL  MCDOFF 

REAL  MOOT 

DIMENSION  ACTI<3/3) 
DIMENSION  CDPOFF  <21 ) 
DIMENSION  CPTC8/2) 


8URN0UT 

LNEXIT 

/ 

EOL 

0 

IMPACT 

0 

LRATE 

SET90 

/ 

SKIP 

IXO 

/ 

I X  F 

0 

IYO 

0 

IYF 

IX 

0 

IY 

0 

II 

IXDOT 

0 

IYDOT 

0 

IZDOT 

0 

ISP 

MCDOFF 

MOOT 

0 

MCDON 

0 

MCMQ 

0 

MCNA 

NOTIME 


DIMENSION  CDPOFF  <21 )  /  CDPON<21>  /  CM8TCI7)  /  CNAT (7) 

DIMENSION  CPTC8/2)  /  CPAF<2) 

DIMENSION  CKC4/50) /  IDAY<3> 

DIMENSION  DCM<3/3>  /  DTDHC1 2> 

DIMENSION  HGcOPB  <1 3) 

DIMENSION  LA9ELC8) 

DIMENSION  MCDOFF  <21 )  /  MCDONC21 )  /  MCMQ<18)  /  MCNA<7) 

DIMENSION  MCP<8) 

DIMENSION  PRESBC12) 

DIMENSION  SPINT < 30)  /  SPNRAT<30> 

DIMENSION  TEMPMB <1  2)  /  THRUST<9)  /  TTIMEt 9) 

DIMENSION  WUVW(3/3)  /  WXYZ<3/3> 

DIMENSION  XX( 3/50) 

DIMENSION  Y<3/50) 

SPNRAT  IS  THE  2.75  AVERAGE  ROLL  RATE  IN  RPS  AS  OF  30  MAY  30 
FROM  MEASURED  NAVY  DATA  PREPARED  BY  D. BROWN.  THE  VALUES  IN  SPINT 
ARE  THE  TIME  BREAKPOINTS  FOR  SPNRAT. 

DATA  SPINT/  0.0/  0.1/  0.2/  0.3/  0.4/  0.5/  0.6/  0.7/  0.8/  1.0/  2.0/ 


/  SPNRAT (30) 


/  TTIMEC9) 


B  9.0/  10.0/  12.0/  16.0/  20.0/  25.0/ 

DATA  SPNRAT/  0.0/  17.0/  28.0/  34.0/  39.0/  43.5/  46.3/  48.7/  51.5/ 
A  54.0/  19.0/  14.0/  8.3/  6.7/  8.5/  11.0/  12.0/  12.2/  11.0/  9.2/ 

8  6.0/  3.6/  0.0/  ”1.0/  ”2.0/  ”2.5/  ”2.5/  ”2.6/  “3.0/  ”3.0/ 

POWER  OFF  DRAG  C OE F F IC I ENT-C DPOF F” FUNCTI ON  OF  MACH 
DATA(CDP0FF<I)/I=1/20)/ 

A  0.700/0.700/0.730/0.809/0.863/0.960/0.977/0.939/1 .000/0.1008/ 

B  1 .010/1 .012/1.005/0.990/0.970/0.940/0.875/0.811/0.765/0.730/ 
POWER  OFF  DRAG  COEFFICIENT  MACH  NUMBER  T ABLE-MCDOFF 
DATACMCDCFFCI)/  1=1/  20)/ 

A  0.0/  0.76/  0.82/  0.90/  0.94/  1.0/  1.03/  1.06/  1.10/  1.15/ 

8  1.18/  1.28/1.34/  1.48/  1.58/1.71/1.94/  2.20/  2.40/2.6/ 

POWER  ON  DRAG  COEFFICIENT-CDPON-FUNCTION  OF  MACH 
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OATA<COPON<I)/  1=1x20) / 

A  C. 5 5/ 0.55/0. 576/0. 629/0. 65/. 6 35/0.699/0. 71/0. 727/0. 741 5/0.747/ 

8. 7 6x0. 75 7x0. 75 3x0.742/ .7 24/0. 6 81x0. 6 5x0. 62 8x0. 61 2/ 

POWER  ON  0 RAG  COEFFICIENT  MACH  NO.  TA8LE-MC00N 
OATA<MCOON<I)x  I  =1/  20)/ 

A  0.0/. 76/  0.82/  0.90/. 94/  1.0/  1.03/  1.06/  1.10/  1.15/  1.18/ 

S  1.28x1.34/1.48  /  1.58/  1.71/  1.94/  2.20x2.40x2.60/ 

CLP- ROLL  DAMPING  COEFFICIENT 
DATA  CLP/-29.539/ 

PITCH  DAMPING  COEF  FICIENT-CMQT-PER  RAD /SEC /-FUNCTION  OF-  MACH 
OATACCMOT  <I) /  1=1x17)/ 

A  -230.5/  -230.5/  -261.0/  -290.3/  -301.8/  -308.1/  -313.2/  -340.0/ 

B  -379.4/  -368.0/  -325.9/  -294.1/  -253.4/  -94.2/  -36.7/  -31.3/ 

C  -30.8/ 

PITCH  OAMPING  MACH  NO.  TABLE-MCMO 
OATA<MCMQ<I) /  1=1/  18)/ 

A  0.0/  0.60/  0.70/  0.80/  0.90/  0.95/  1.00/  1.05/  1.10/  1.20/ 

B  1.30/  1.40/  1.50/  2.00/  3.00/  4.0/  5.0/  10.0/ 

NORMAL  FORCE  COEFFICIENT  SLOPE-PER  D EG  RE  E-CNAT-F UNCTION  OF  MACH 
OATACCNATCI)/  1=1/  7)/ 

A  12.44x12.44x12.597x8.874x6.585x6.012/6.012/ 

NORMAL  FORCE  MACH  NO.  TA3LE-MCNA 
D ATA (MCN A ( I ) /  1=1/  7)/ 

A  0.0/  0.60/  0.80/  1.75/  2.50/  3.00/  10.0/ 

ANGLE  OF  ATTACH  TA9LE  FOR  CP-IN  OEGREES-CPAF 
DATACCPAFCI)/  1=1/  2)/ 

A  0.0/  90.0/ 

CENTER  OF  PRES SURE-CPT-F UNCTION  OF  MACH  AND  ANGLE  OF  ATTACK 

IN  CALIBERS  FROM  NOSE 

OATACCPT Cl/1 )/  I  =  1x7)/ 

A  14.54/  14.54/  14.75/  15.29/  15.07/  13.31/  12.21/ 

B  <CPT<1/2)x  1=1x7)/ 

C  14.54/  14.54/  14.75/  15.29/  15.07/  13.31/  12.21/ 

CENTER  OF  PRESSURE  MACH  NO.  TABLE-MCP 
OATACMCPCI)/  1*1/  8)/ 

A  0.0/  0.60/  0.80/  1.00/  1.75/  2.50/  3.00/  10.0/ 

INITIAL  AND  FINAL  VALUES  OF  CENTER  OF  GRAVITVCGO/CGF-FEET  FROM  NOSE 
DATA  CGO/  CGF/2.11 6667/  1.800000/ 

REFERENCE  01 AMETER-0 REF-2.75  INCHES-IN  FEET 
DATA  OREF/0. 2292/ 

ATMOSPHERIC  OENSITY  REFERENCE  TABLE 
OATACDTOHU)/  1=2/12)/ 

A  -0.0065/  0.0/  0.003/  0.0/  0.004/  -0.0045/  0.0/  0.02/  0.31/ 

C  0.005/  0.0035/ 

DATA  F/0. 34838395/ 

OATA<HGEOPB<I)x  1=  2/  13)/ 

A  0.0/  11000./  25000.x  47000./  53000./  79000./  90000./  105000.x 
8  160000.x  170000.x  200000.x  700000./ 

1  SLU6=  32.1739  LB.  OF  MASS  (L9-S9. IN. =0.00021 58  SLUG-SQ.FT.) 

INITIAL  ANO  FINAL  AXIAL  MOMENT  OF  INERTI A-IXO/ I XF- SLUG-F EET*F EET 
OATA  1X0/  IXF/5 .6755E-3/  2.4621E-3/ 

INITIAL  ANO  FINAL  TRANSVERSE  MOMENT  OF  INERTIA-I YO/ IYF 
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DATA  IYO/IYF/1. 4449968*  1.1869/ 


DATA  IMPACT/. FALSE./ 

DATA  NOTIME/. FALSE./ 

DATA<PRESB<I)/  I  =  2*12)/ 

A  1 .0332E4/  2.3078E3/  2.5381  E 2 1.2285E1/  5.9478/  1.029E-1/ 
B  1.066E-2/  7.618E-4/  3.691E-5/  2.88E-5/  1.454E-5/ 

DATA  PI/3.14159/ 

DATA  RADCON/57. 2957795/ 

C  REFERENCE  SURFACE  AREA  OF  EARTH 
DATA  S/6.378178E+6/ 


C  MOTOR  BURN  TINE-SUSOFF 
DATA  SUSOFF/1 .075/ 

C  ATMOSPHERIC  TEMPERATURE  REFERENCE  TABLE 
DATA<TEMPM8<I)/  1*2/  12)/ 

A  288.16/  216.66/  216.66/  282.66/  282.66/  165.66/ 

B  165.66/  225.66/  1325.66/  1425.66/  1575.66/ 

c  thrust  table-thrust-pounds;time  TABLE  OF  THRUST 

DATACTHRUSTCI)/  1=1/  9)/ 

A  0.0/  1600.0/  1350.0/  1300.0/  1550.0/  1750.0/  1550.0/  0.0/  0.0/ 
DATACTTIME(I)/  1=1/  9)/ 

A  0.0/  0.05/  0.2/  0.3/  0.8/  0.95/  1.00/  1.075/  10.0/ 

C  TOTAL  IMPULSE-TI 
DATA  TI/1510./ 

C  INITIAL  AND  FINAL  VALUE  OF  MISSILE  WEIGHT  WO/  WF 
DATA  WO/  WF/23.75/  16.59/ 

C  WEIGHT  WHEN  ROUND  CLEARS  LAUNCHER-WL 
DATA  WL/23.296/ 

C  MISSILE  LENGTH  IN  FEET-XLENG 
DATA  XLENG/5 .491 67/ 

C  FIND  THE  CURRENT  DATE  AS  DATE (IDAY) 

CALL  DATECID AY) 


DATA  STATEMENTS  INITIALIZING  ALL  VARIABLES  IN  NAMELIST  PUTT 
DATA  ALTIN/-50.0/ 

DATA  DELTT/  DEL02/  DIVANG/O. 001 /  0.0005/  0.0/ 

DATA  DNSCH/  DPIT/  DT/  DYAW/  DZZ/0.5/  0.0/  0.001/  2*0.0/ 
DATA  EOL/. FALSE./ 

DATA  FACTOR/  FC A/  FTORM/  GREFF/0.0/  -6.56/  0.0/  32.174/ 
DATA  HELIC/  N/  NN/O.O/  15/  1/ 

DATA  PHI  ZD/  PITHAL/  PSIZD/  QELNCH/  3*0.0/  16.0/ 

DATA  RAMP/  REF/  RK/  ROTRAD/O.O/  20902703.0/  4.0/  22.0/ 
DATA  SCARF/0.0/ 

DATA  TDC/THTZD/TFINAL/TIME/3.0/0./7.0/.09/ 

DATA  TLL/TPRINT/TROTZ.09/.1/.092/ 

DATA  VINIT/  WCF/  WDF/68.0/  0.0/  0.0/ 

DATA  XMV/  XT AR/1 46.5/  13000.0/ 

DATA  YAWMAL/  YTAR/O.O/  0.0/ 

DATA  ZTAR/O.O/ 

WRITE (6/ 3000)  IDAY  54 
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WRITE  <6,  36)  ALTIN,DELTT,DEL02,DIVANG 
WRITE<6,  37)  DNSCH,DPIT,DT,DYAW 
WRIT  E(6, 38)  DZZ,£OL, FACTOR, FCA 
WRITE<6,39)  FTORM, GREFF, HELIC, N 
WRITE  <6,  40)  NN,PHIZD,PITMAL,PSIZD 
WRITE(6,41)  QELNCH, RAMP, REF, RK 
WRITE (6x42)  ROTRAD,SCA«F,TDC,THTZD 
WRITE <6, 43)  TFINAL,TIM£,TLL,TPRINT 
WR ITE  < 6, 44)  TROT, VINIT, WCF,WDF 
WRITE(6,45)  XMV,XTAR,YAWMAL,YTAR 
WRITE ( 6, 46)  ZTAR 

vinit  is  the  airspeed;xtar  IS  THE  range;qelnch  IS  THE  LAUNCHER 
ELEVATION  (ANGLE  BETWEEN  LAUNCHER  CENTERLINE  AND  HELICOPTER 
CENTERLINE>;XHV  IS  THE  MISSLE  VELOCITY  AT  LAUNCHER  EXIT;  TIME, 

TLL  ARE  VARIABLES  OF  TIME  AT  LAUNCHER  EXIT,  HELIC  IS  THE  DOWNWASH 
INDICAT0R:1-D0WNWASH  USED,0=N0  DOWNWASH  USEDjFTORM  IS  UNIT  USAGE: 
ENGLISH  SYSTEM  =0, METRIC  SYSTEM=1. 

NAMELIST/PUTT/ 

A  ALTIN, 

B  DELTT,  DEL02,  DIVANG, 

C  DNSCH,  DPIT,  DT,  DYAW,  DZZ, 

D  EOL, 

E  FACTOR,  FCA,  FTORM,  GREFF, 

F  HELIC,  N,  NN, 

H  PHI  ZD,  PITMAL,  PSIZD,  QELNCH, 

I  RAMP,  REF,  RK,  ROTRAD, 

J  SCARF, 

K  TDC,  THTZD,  TFINAL,  TIME, 

L  TLL,  TPRINT ,  TROT, 

M  VINIT,  WCF,  WDF, 

N  XMV,  XTAR, 

0  YAWMAL,  YTAR, 

P  ZTAR 

READ ( 5, PUTT) 

C  THE  DERIVATIVES  ARE  FOUND  IN  THE  XX(1,J)  SECTION  OF  THE  XX  ARRAY. 

C  XX(2, J),S  ARE  INTEGRALS  OF  XX(1,J),S.  THE  FOLLOWING  IS  A  LIST  OF 
C  THE  DERIVATIVES  8Y  XX(1,J)  LOCATIONS. 

C  XX(1,1)-ANGULAR  ACCELERATION  ABOUT  THE  MISSILE  X-AXIS. 

C  XX<1,2)-ANGULAR  ACCELERATION  ABOUT  THE  MISSILE  Y-AXIS. 

C  XX(1,3)-ANGULAR  ACCELERATION  ABOUT  THE  MISSILE  Z'AXIS. 

C  XX <1 ,4) -LINEAR  ACCELERATION  ALONG  THE  MISSILE  X-AXIS. 

C  XX(1,5)-LINEAR  ACCELERATION  ALONG  THE  MISSILE  Y-AXIS. 

C  XX(1, 6) -LINEAR  ACCELERATION  ALONG  THE  MISSILE  Z-AXIS. 

C  XX<1 ,7)-LINEAR  VELOCITY  ALONG  THE  GROUND  X(RANGE)-AXIS. 

C  XX(1,8)-LINEAR  VELOCITY  ALONG  THE  GROUND  Y<DEFLECTION)-AXIS. 

C  XX(1,9)-LINEAR  VELOCITY  ALONG  THE  GROUND  Z(VERTICAL)-AXIS. 

C  XX(1,10)-DERIVATIVE  OF  THE  EULER  ANGLE  SI  =  PSI 

C  XXd, ID-DERIVATIVE  OF  THE  EULER  ANGLE  TH  =  THETA 

C  XX(1, 12) -DERIVATIVE  OF  THE  EULER  ANGLE  FI  =  PHI. 

C  XX(1,13)-DERIVATIVE  OF  THE  MASS. 

C  XX(1,14)-ACCELERATI0N  ALONG  LAUNCHER  CENTERLINE. 

C  XXd, 15) -  VELOCITY  ALONG  LAUNCHER  CENTERLINE. 

C  THIS  IS  THE  INITIAL  SECTION  OF  THE  PROGRAM  AND  IS 

C  ONLY  PASSED  THRU  AT  THE  BEGINNING  OF  EACH  RUN. 

C  SET  ARRAYS  TO  ZERO 
00  50  I  *  1,  3 
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C  INITIALIZE  BURNOUT  VARIABLE 

BURNOUT  =  .FALSE. 

C  WHEN  THE  VARIABLE  LRATE  IS  TRUE, THE  ROUND  CG  HAS  CLEARED  THE 

C  LAUNCH  TUBE. AT  THIS  TIME, THE  MALLAUNCH  RATES  PITMAL  AND  YAWMAL  ARE 

C  ADDED  TO  THE  BODY  RATES. ASSUMPTION-ROUND  LENGTH  XLENG= LAUNCHER 
C  LEN6TH. 

LRATE  =  .FALSE. 

C  INITIALIZE  VARIABLE  THAT  SETS  MISSILE  PARAMETERS  TO  THER  VALUES 

C  WHEN  BURNOUT  OCCURS. 

SET80  *  .FALSE. 

C  INITIALIZE  L06ICAL  VARIABLE  THAT  ALLOWS  PROGRAM  TO  SKIP  SETTING 

C  VALUES  AT  BURNOUT  WHEN  THE  LOGICAL  SETBO  IS  TRUE. 

SKIP  =  .FALSE. 

C  CALCULATE  TIME  RATE  OF  CHANGE  OF  CENTER  OF  GRAVITY-CGDOT 

CGDOT  =  (CGO  -CGF) /SUSOF F 
C  INITIALIZE  CENTER  OF  GRAVITY-CG 

C 6  *  CGO 

C  CALCULATE  MASS  AT  TIME  ZERO 

XMASS-WO/GREFF 
XX(2,1 3)  =  XMASS 

C  CALCULATE  TIME  RATE  OF  CHANGE  IN  AXIAL  AND  TRANSVERSE  MOMENTS 

C  OF  INERTIA-IXDOT,  IYDOT 

IXDOT  =  <1X0  -  IXFJ/SUSOFF 
IYDOT  =  (I YO  -  IYF)/ SUSOFF 
IZDOT  =  IYDOT 

C  CALCULATE  SPECIFIC  IMPULSE-1 SP=T0TAL  IMPULSE/PROPELLANT  WEIGHT 

ISP  =  TI / ( WO  ~  WF) 

C  SET  FLAG  FOR  SUBROUTINE  HELVEL 

IFL  =  0 

C  CALCULATE  QUADRANT  ELEVATION  IN  MILS  AND  AIRSPEED  IN  KNOTS 

QEMILS  =  QELNCH*10G0./RA0C0N 
AIRSP  *  VINIT*3600./60S0. 

C  CALCULATE  HELICOPTER  ATTITUDE-HELAT-IN  DEGS 

HELAT  *  0.0 

IF<AIRSP  .GE.  50.)  HELAT  =  FCA*(AIRSP  -  50.J/70. 

C  CALCULATE  REFERENCE  AREA-AREFF 

AREFF  =  <PI*DREF*DREF)/A. 

C  SET  MISSILE  ALTITUDE  XX < 2, 9) =INITI AL  ALTITUDE  ALTIN 

XX<2,9)  =  ALTIN 

QELNCH  -  ANGLE  BETWEEN  LAUNCHER  CENTERLINE  AND  HELICOPTER 
CENTERLINE 

HELAT  -  ANGLE  BETWEEN  HELICOPTER  CENTERLINE  AND  LOCAL  HORIZONTAL 
THL  -  ANGLE  BETWEEN  LAUNCHER  CENTERLINE  AND  LOCAL  HORIZONTAL 
DIVANG  -  ANGLE  BETWEEN  HELICOPTER  VELOCITY  VECTOR  VINIT  AND  LOCAL 
HORIZONTAL 

ALV  -  ANGLE  BETWEEN  LAUNCHER  CENTERLINE  AND  HELICOPTER  VELOCITY 
VECTOR  VINIT 

THL=  (HELAT  *  QELNCH)/RADCON 
ALV  =  THL“<DIVANG/RADC0N) 

C  CALCULATE  VELOCITIES  IN  THE  BODY  X  AND  Z  DIRECTION  DUE  TO 

C  HELICOPTER  VELOCITY  VINIT. 

XX<2,4)  =  VI NI T*  CO  S ( ALV) 

XX(2,6)  =  VINIT*  3IN<  ALV) 

C  SET  THE  EULER  ANGLES  TO  THEIR  INITIAL  VALUES.  XX<2,10>  =  SI, 

C  XX<2,11)=TH,XX<2,1 2)=FI. 

XX(2,10)  =  PSIZD/RADCON 

XX  <2,1 1 )  =  THL 

XX(2,1 2)  *  PHIZD/RADCON 

C  IF  EOL (END  OF  LAUNCH)VARIA3LE  IS  TRUE,TR AJ ECTORY  BEGINS  WITH 

C  PHYSICAL  PARAMETERS  AT  THE  TIME  THE  ROUND  CLEARS  THE  TUBE. 

I F < . NOT ■  EOL)  GO  TO  60 
C  XMV*MISSILE  VELOCITY  AT  LAUNCH 


56 


nonnnnnnono 


XX(2/4)  =  VINIT*COS(ALV)  ♦  XMV 
C  WL-WEI6HT  OF  ROUND  AT  TU3E  EXIT. 

XMASS  =  WL/6REFF 
XX(2s1 3)  =  XMASS 

C  SET  THE  BODY  RATES  Q  AND  R  TO  THE  MALLAUNCH  RATES 
XX(2s2)  =  PITMAL 
XX (2/3)  =  YAWMAL 

C  SET  LRATE  TRUE  TO  INOICATE  MALLAUNCH  RATES  HAVE  BEEN  ADDED 
LRATE  =  .TRUE. 

C  CALCULATE  CENTER  OF  GRAVITY  AT  END  OF  LAUNCH 

CG  =  C  GO  “  TIME*C6D0T 
60  CONTINUE 

C  INITIALIZE  OTHER  PARAMETERS 
IPRNT  =  1 

TTO  =  SUSOFF  +  DEL02 
TEXIT  =  TLL  +  DEL02 

C  DXLN/  D YLN/  DZLN  REPRESENT  THE  XYZ  COORDINATES  OF  THE  CG  OF  THE 

C  LAUNCHER  W/R  TO  THE  HELICOPTER  ROTOR  HUB. 

DXLN  =  2. 

D YLN  =  3.5 
DZLN  =  7. 

DXX  =  0. 

DYY  =  0. 

END  OF  INITIALIZATION  SECTION  OF  MAIN  PROGRAM 

READ  ONE  CARD  TO  DESCRIBE  PERTINENT  FEATURES  ABOUT  PARTICULAR 
TRAJECTORY  RUN. 

30  FORMAT (8 Al 0) 

CHANGED  READ/WRITE  LABEL  TO  COMMENT 
READ(5/30)  LABEL 

31  FORMAT (1  HI ) 

WRITE ( 6/ 31 ) 

UR  IT  E ( 6s  30)  LABEL 

32  FORMAT (1  HO) 

3000  F0RMAT(60XsI2s,/'sI2s'/,sI2s) 

WRITE(6s32) 

WRITE(6s25) 

25  FORMATt  ***********  2.75  MISSLE  WITH  Ml  51  WARHEAD*) 

34  F0RMAT(20Xs' TRAJECTORY  RUN  WAS  MADE  WITH  HELICOPTER  DOWNWASH  PROFI 
1LE' ) 

35  FORM AT ( 20Xs '  $$$  NO  HELICOPTER  DOWNWASH  PROFILE  WAS  USED  IN  THIS  T 
1  RAJ ECT OR Y  RUN*) 

IF (HELIC  .E Q.  1.)  WRlTE(6s34) 

IFCHELIC  . E 0 .  0.)  WRITE<6s35) 

C  WRITE(6sPUTT) 

C  THE  GROUP  OF  WRITE  STATEMENTS  THAT  FOLLOW  CAUSE  A  PRINTOUT  DOWN 
C  THE  PAGEsIF  ACROSS  THE  PAGE  GO  TO  THE  PRECEEDING  WRITE/ "WRITE (6s PUTT) 
WRITE (6s 36)  ALTINs  DELTTs  DEL02sDI V ANG 

36  FORMAT (  *  ALTIN=' s F6. 2s2Xs ' DELTT*1 s E«. 4s 2Xs * DEL02*' s E9.4 s2Xs *  DIVAN 
1G='sF5.2) 

WRITE(6s37)  DNSCHsDPITsDTsDYAW 

37  FORMAT (  *  DNSCH='sF5.2s2<s,DPIT=*sF5.2s2Xs'DT=*sE9.4s2Xs'DYAW=*s 
1F5.2) 

WRITE(6s38)  DZZsEOLsFACTORsFCA 

38  FORMAT {  •  DZ Z  =  ' s F5 . 2s 2Xs * EOL  =  * sL3s 2Xs *  FACTOR** s F5. 2s 2Xs • FC A=' s 
1F8.3) 

WRITE ( 6s  39)  FTORMsGREFFsHELICsN 

39  FORMATt  *  FTORM*’ sF5 . 2s2 Xs * GRE FF=* s F 8. 3s2Xs • HELIC*' s F3 . 2/2Xs ' N*' s 
1110.4) 

WRITE<6s40)  NNsPHIZDsPITMALsPSIZD 

40  FORMAT <  *  NNa'sI10.2s2Xs,PHIZD  =  ,sF5.2s2Xs'PITMAL  =  ,sF5.2s2Xs'PSIZD* 
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1*,F5.2) 

WRITEC6/41)  QELNCH, RAMP, REF, RK 

41  FORMAT <  '  9£LNCH='/F6.1/2X,'RAMP=',F5.2/2X/'REF='/E12.4/2X/'RK='/ 
1F5.2) 

WRITEC6/42)  ROTRAD/SCARF/TDC/THTZD 

42  FORMAT  <  •  ROTR  AD  =  *  ,  F6. 2/ 2X  ,  *  S  C  AR  F=  •  ,  F  5 . 2,  2X,  '  TOC  =  '  ,  F  5  .  2,  2  X,  •  THTZD  = 
1 '/F5.2) 

WRITE <6/ 43)  TFINAL/TIME/TLL/TPRINT 

43  FORMAT <  '  TFINAL=*,F5.2/2X/*TIME='/E9.4,2X/'TLL='/£9.4/2X/'TPRINT= 
1',F5.4) 

WRITEC6/44)  TROT/VINIT/WCF/WDF 

44  FORMAT (  •  TROT='/E?.4/2X/'VINIT='/F6.2/2X/*WCF='/F5.2/2X/'WDF='/ 
1F5.2) 

WRITEC6/45)  XMV/XTAR/YAWMAL/YTAR 

45  FORMATC  *  XMV='/F7.2/2X/'XTAR=',F9.2/2X/*YAWMAL='/F5.2/2X/'YTAR='/ 
1F5.2) 

WRITEC6/46)  ZTAR 

46  FORMAT  <  '  ZTAR=*F5.2) 

WRITE <6/  33)  DIVANG,  HELAT/  QEMILS,  AIRSP,  WO 
33  FORMATC  SIX  DEGREE  OF  FREEOOM  SIMULATION  FOR  2.75  ROCKET  WITH  H 
1 ELICOPTER  DOWNWASH  PROFILE'//*  HELICOPTER  DIVE  ANGLE  =  '/F7.2,' 
2DEGREES* /•  HELICOPTER  ATTITUDE  =',F6.2,'  DEGREES'/'  LAUNCHER  9 
3.E.  =  ' / F8. 2/ *  MILS'/'  INDICATED  AIRSPEED  =  '/F7.2,'  KNOTS'/' 
4INITIAL  ROCKET  W EIGNT= '/ F6 . 2, ' POUNDS '/ ) 

C  THE  PROGRAM  NOW  OFFERS  A  CHOICE  8ETWEEN  2ND  ORDER  AND  4TH 

C  ORDER  RUNGE  -KUTTA  INTEGRATION.  IF  RK=2.0,2ND  ORDER  RUNGE-KUTTA 
C  INTEGRATION  IS  USED.  THE  DEFAULT  VALUE  OF  RK  IS  4.0*  AND  CAN  3E 
C  CHANGED  THRU  NAMELIST  PUTT. 

K2  =  2 
<=4 

T0=TIME 

85  IF<K2  -  2)2/9,2060 
7  IFCK-4)2/9,2060 

2060  WRITEC6, 2061 ) 

2061  FORMATC  ERROR  -K  EXCEEDS  4') 

9  SI  =  XX( 2/10) 

TH  =  XXC2/11) 

FI  =  XXC2/12) 

SSI  =  SIN(SI) 

CSI  =  COS(SI) 

STH  =  SIN(TH) 

CTH  =  COSCTH) 

SFI  *  SIN(FI) 

CFI  *  COSCFI) 

C  CALCULATE  THE  OCM  MATRIX/THAT  TRANSFORMS  FROM  GROUND  TO  BODY/ 

C  USING  THE  CURRENT  VALUES  OF  THE  EULER  ANGLES. 

DCMC1/1)  =  CTH*CSI 

DCMC1/2)  =  CTH*SSI 

DCM(1/3)  =  -STH 

0CM<2/1 )  =  STH*CSI*SFI  -  SSI*CFI 

DCM<2/ 2)  *  STH*SSI*SFI  ♦  CSI*CFI 

DCNC2/3)  *  CTH*SFI 

DCNC3/1 )  =  STH*CSI*CFI  ♦  SSI*SFI 

DCM<3/2)  =  STH*SSI*CFI  -  CSI*SFI 

DCM( 3/ 3)  *  CTH*CFI 

C  CHANGE  INTEGRATION  DELTA  T  AT  3.0  SECONDS 

IFCTIME  .GE.  TDC)  DELTT  =  DT 
DEL02  »  0.5*DELTT 

C  COMPUTE  THE  COMPONENTS  OF  THE  GRAVITATIONAL  ACCELERATION  ALONG 

C  THE  MISSILE  XYZ  AXES  A  G AX/  GAY/  GAZ. 

GAX  =  DCMC1,3)*GREFF 

6AY  *  DCM<2/3)*GREFF  58 


GAZ  =  DC  M ( 3x  3) *G  RE  FF 
C  ORTHOGONALITY  CONSTRAINT 

CALL  INVERT(DCMxACTI) 

DO  84  1=1x3 
DO  84  J=1 x3 

34  DCM<IxJ)=<ACTI(JxI)+DCM<IxJ))/2. 

2  IF<RAMP.NE.O.)WCF=RAMP*TIME 

C  CALCULATE  ROLL  RATE* IN  REV/SECxAS  FUNCTION  OF  TIME  AS  THE  VARIABLE 

C  ROLL. 

ROLL  =  T ABLK (SPNRATx  SPINTx  TIMEx  30) 

XX(2x1 )  =  R0lL*2.*PI 
IF(.NOT.LNEXIT)  XX<2x1)  =  0.0 

C  THE  COMPONENT  OF  THE  TOTAL  ROUND  VELOCITY  VECTOR  ALONG  THE  MISSILE 

C  X  AXIS  IS  U=XX (2x4 ) .  U  IS  FOUND  BY  INTEGRATING  ALL  THE  LINEAR 

C  ACCELERATION  TERMS.  WHAT  IS  DESIRED  IS  ONLY  THE  LINEAR  ACCELERATIONS 
C  ALONG  THE  X  AXIS  THAT  ARE  DUE  TO  THRUSTx  DRAGx  AND  GRAVITYx  IN 
C  ORDER  THAT  WE  CAN  INTEGRATE  TWICE  AND  OBTAIN  THE  DISTANCE  THE  ROUND 
C  HAS  TRAVELED  ALONG  THE  LAUNCHER  CENTERLINE.  LET  UDOTL  3E  THE  ACCEL. 

C  ALONG  THE  MISSILE  X  AXIS  DUE  ONLY  TO  THRUSTxDRAGxAND  GRAVITY. 

UDOTL  =  UMF1  ♦  G AX 

C  ESTABLISH  UDOTL  AS  A  DERIVATIVE. 

XX (1 xl 4)  =  UDOTL 

C  THE  STATExOR  VELOC ITYx OBTAINED  BY  INTEGRATING  XX(1x14)  IS  XX<2x14) 

C  ESTABLISH  THE  DERIVATIVE  OF  DISTANCE  ALONG  THE  LAUNCHER  CENTERLINE 
C  THE  VELOCITYx  AS  XXCIxIS). 

UL  =  XX (2x14) 

XX(1x15)  =  UL 

C  THE  DISTANCE  THE  ROUND  HAS  TRAVELED  ALONG  THE  LAUNCHER  CENTERLINE 

C  IS  THE  STATE  XX(2x15)xTHE  INTEGRAL  OF  XX(1x15>. 

CGTRVL  =  XX( 2x1 5) 

C  THE  ROUND  WILL  HAVE  CLEARED  THE  LAUNCHER  WHEN  THE  CENTER  OF 

C  GRAVITY  HAS  MOVED  A  DISTANCE  =  THE  LENGTH  OF  THE  ROUNDx XLENG .  THIS 

C  WILL  BE  SIGNALED  BY  THE  LOGICAL  VARIABLE  LNEXIT  CHANGING  FROM  FALSE 

C  TO  TRUE. 

LNEXIT  =  CGTRVL  .GE.  XLENG 

C  THIS  SECTION  ADDS  THE  MALLAUNCH  RATE  WHEN  THE  ROUND  CG  HAS 

C  CLEARED  THE  LAUNCHER. ALSOx ADDS  THE  HELICOPTER  VELOCITY  IN  X  AND  Z 
C  AXES. 

IF(LRATE)  GO  TO  2062 
I F (  .NOT. LNEXIT)  GO  TO  2062 
LR ATE  =  .TRUE. 

XX(2x2)  =  XX(2x2)  ♦  PITHAL 
XX (2x  3 )  =  XX ( 2x  3 )  +  YAWMAL 
JUMP=3 
GO  TO  8 

2062  REZIF=REF-XX<2x9) 

RF=SQRT(XX(2x7)**2+XX(2x8)**2+REZIF**2) 

HEIGHT=RF-REF 

IF(HEIGHT)4xSx5 

4  IF(HEIGHT+1.)3x6x6 
6  HEIGHT=0. 

GO  TO  5 

3  WRITE<6x1 501 ) 

JUHP=2 

GO  TO  8 

5  SQIRT=SQRT<XX<2x7)**2+XX<2x8)**2) 

TANG  =  SQI RT/R  F 

RGF*REF*ATAN(TANG/SQRT(1.-TANG**2)) 

WIND=SQRT<WCF*WCF+WDF*WDF> 

IF(. NOT. LNEXIT)  GO  TO  230 
IFCHELIC.EQ.O.)  GO  TO  230  59 

DXX=Oa 


DZZ»0. 

I F (TIME. 6T . . 3)  SO  TO  230 
TPHASE=TIHE+TROT 
XH=VINIT*COS  (DIVANG)*TIHE+DXLN 
YH=DYLN 

ZN*VINIT*SIN(DIVANG)*TIME-ALTIN+DZLN 

XHXH*XV(2,7>-XH 

YMYH=kX(2,8)-YH 

ZHZH=XX(2,9)+ZH 

DDIS=SQRT(XHXH**2+YMYH**2+ZMZH**2> 

THAT*DI VAN6+HELAT 
THEO-O. 

I F ( ZHZH. NE . 0 • )  THEO=ATAN(XNXH/ZMZH> 

THTH=THEO-TH AT 
XD=-DDI$*SIN(THTH) 

ZD=-DDI$*COS(THTH> 

YD=-YH 

DDI$S=$QRT(XD*XD+YD*YD> 

IF(DDISS.LE. ROTRAO) CALL  HELV£L(XD, YD/ZD,DWX/DWY/DWZ,TPHASE/I FL) 
I  FL=1 

DXX=DWX*COS(THAT)-DWZ*SIN(THAT> 

OYY=DWY 

DZZ=“DHZ*COS (THAT) “DyX*S IN (THAT) 

230  CONTINUE 

IF(SQIRT.6T. ROTRAO)  DZZ=0. 

IF(WIND.GT.O..AND.SQIRT.GT.O.)  GO  TO  10 
WXYZd/1  )=WDF+DXX 
WXYZ(2,1)=WCF+DYY 
HXYZ(3,1)=DZZ 
GO  TO  12 
10  CONTINUE 

yXYZ(1,1)=(-yCF*XX(2,8)+yDF*XX(2,7)*REZIF/RF>/SQIRT+DXX 
WXYZ(2,1)=(yCF*XX(2,7)+WDF*XX(Z,3>*REZIF/RF)/S9IRT+DYY 
WXYZ(3,1)=MDF*SQIRT/RF+DZZ 
12  CALL  MULTY(DCH,WXYZ,1,yUVW) 

IF(FTORM)223,222,223 

222  HEIGHT=HEIGHT/3. 2808333 

223  HGTIND=0.0 
H6E0P=S*HEIGHT/(S+HEIGHT> 

I F(HGE0P-90000. >201,201, 202 

201  UM0L=28. 966 
GO  TO  205 

202  IF(HGEOP-1SOOOO. >203,203/204 

203  WHOL=22.0-5.04483574*ATAN(6356.766*.04*HEIGHT/(6356.766+ 
1HEI6HT>-8.8> 

GO  TO  205 

204  WM0L=27. 1 06-7. 9356971 *ATAN((HGE0P-1 80000.0) /1 40000.0) 

205  1=1 

206  1=1+1 

I  F(HGEOPBd)-HGEOP)  207,209, 210 

207  IF(I-13)206,208,208 

208  HGTIND=1.0 
GO  TO  15 

209  TEMPN0=TEMPN3(I) 

H6E0P0=HGE0P8(I) 

OOOH=OTOH(I) 

PRES0*PRES8( I) 

SO  TO  211 

210  000H*0T0H(I“1) 

PRES0=PRESB(I-1)  60 

H6EOPO*HSEOP3( I”1 ) 


TEMPM0=T£MPMB(I-1) 

21 1  TEMPM=TEMPMO+DODH*(HGEOP-HGEOPO> 

TEMP =TEMPM*WM0L/28. 966 
IF(HEIGHT-9C 000. >216/216/217 

216  VS0UND=20.04633*SGRT(TEMPM) 

60  TO  218 

217  VS0UNDs20.04633*SQRT(165.66) 

218  IF(00DH)212/213/212 

212  PRES=PRESO*(TEMPMO/(TEMPMO+DODH*(HGEOP-HGEOPO> >>**(. 0341  65/DODH) 
GO  TO  214 

21 3  PRES*PRESO*£XP(-.034165*(HGEOP-HGEOPO)/TEMPMO> 

214  RH0=(F*PRES)/(TEMPM*10.1971> 

I F (FTORM .NE. 0. )  GO  TO  15 
V$0UNDs3.2303333*VS0UND 

NOTE  DENSITY  IS  COMPUTED  IN  NEWT0NS/M**3  OR  SLU6S/FT**3 
RHO=RHO/5l5.375 
HEIGHTS.  28033  33  *H  EIGHT 
15  IF(HGTIND.LE.O.)  GO  TO  20 
WRITE(6/1 500) 

JUMP=2 
GO  TO  3 

1500  F0RMAT(38H  HEIGHT  EXCEEDS  MAXIMUM  VALUE  IN  TABLE) 

1501  FORMATO  9H  HEIGHT  IS  NEGATIVE) 

20  XMW=XX(2/4)-WUVW(1/1> 

YMW*XX(2/5)-WUVW(2/1 ) 

ZMW=XX(2/6)-WUVW(3/1) 

VVW=YMW*YMW+ZMW*ZMW 

VRF=SQRT(VVW*XMW*XMW> 

XMACHN=VRF/V SOUND 
VALF=S9RT(VVW/(XMW**2)) 

ALFTOT*ATAN(VALF)*RADCON 
IF(.NOT.LNEXIT)  ALFTOT  *  0.0 
BANK*1 .5707 
IF(YMW.EQ.0.)8ANK=0. 

IFCZMW.NE.O. ) BANK=ATAN2(YMW/ZMW> 

IF(XMW.£9.0.)  GO  TO  503 
ALPPRN=ATAN2(ZMW/XMW) 

3ETPRN*ATAN2(YMW/XMW) 

60  TO  510 
503  ALPPRN=0. 

BETPRN*0. 

510  IF (VRF  .NE.  0.)  GO  TO  512 

511  FORMATd OX/' RELATIVE  VELOCITY  VRF  WAS  ZERO-RUN  STOPPED  TO  PREVENT 
1DIVISI0N  BY  ZERO') 

WRITEC6/511) 

STOP 

512  CONTINUE 

CALCULATE  NORMAL  AND  SIDE  FORCE  COEFFICIENTS  CN  AND  CY  AS 
FUNCTIONS  OF  MACH  NUMBER. 

CN  *  TABLK(CNAT/  MCNA/  XMACHN/  7)*ALFT0T 
CY  =  CN 

MODIFY  CN  AND  CY  BY  THE  BANK  ANGLE 
CN  *  -CN*COS (BANK) 

CY  =  -CY*SIN(8ANK) 

CHANGE  BANK  ANGLE  TO  DEGREES  FOR  PRINTING 
BANK  =  8ANK*RADC0N 

CALCULATE  THE  AXIAL  FORCE  COEFFICIENT-C A-AS  A  FUNCTION  OF  MACH 
NUMBER  AND  WETHER  OR  NOT  THE  MOTOR  IS  BURNING. 

IF(BURNOUT)  GO  TO  513 

POWER  ON  DRAG  COEFFICIENT  -  CA 

CA  =  TABLK(CDPON/  MCDON/  XMACHN/  21) 

GO  TO  514 
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513  CONTINUE 

C  POWER  OFF  DRAG  COEFFICIENT  -  CA 

CA  =  T  A8LK ( C DPOF F*  HCDOFF*  XMACHN*  21) 

•liA  CONTINUE 

C  CALCULATE  TH RUST-THRSTP- AS  A  FUNCTION  OF  TIME 

THRSTP  =  TASLtUTHRUST*  TTIME*  TIME*  9) 

C  CALCULATE  THE  AERODYNAMIC  FORCES  ALONG  THE  MISSILE  XYZ  AXES  AS 
C  FUl AP*FV1AP*FW1AP 
C  TP  =  DYNAMIC  PRESSURE 
DP  =DNSCH*RHO*VRF**2 
DPS  =  DP*ARE  FF 
FUlAP  =  -CA*DPS 

FV1AP  =  CY*DPS 

FW1AP  =  CN*DPS 

C  IN  ORDER  TO  ALIGN  THE  80DY  WITH  THE  THRUST*ROTATE  FROM  BODY  TO 
C  THRUST  AXIS  FIRST  THRU  DYAW  ABOUT  Z  AXIS*THEN  THRU  DPIT  ABOUT  NEW 
C  Y  AXIS. 

C  CALCULATE  THE  THRUST  FORCES  ALONG  THE  MISSILE  XYZ  AXES  AS  FU1GP* 

C  FV1GP*  FW1GP. 

FU1GP  =  THRSTP*COS<DPIT)*COS<DYAW> 

FV1GP  *  THRSTP*COS(DPIT)*SIN(DYAW) 

FW1GP  =  THRSTP*SIN(DPIT> 

C  CALCULATE  THE  LINEAR  ACCELERATIONS  ALONG  MISSILE  XYZ  AXES  AS  UMF1 * 

C  VMF1 *  WMF1  DUE  TO  THRUST  AND  DRAG. 

XMASS  *  XX(2*1 3) 

UMF1  *  (FU1AP  +  FU1 6P) /XMASS 
VMF1  «  CFV1AP  ♦  F VI GP) /XMASS 
WMF1  *  CFW1AP  +  FW1GP)/XMASS 

C  DEFINE  THE  LINEAR  AND  ANGULAR  VELOCITIES  ALONG  MISSILE  XYZ  AXES  AS 
C  Us  Vs  Ws  Ps  9s  R. 

U  *  XX (2s A) 

V  *  XX(2s5) 

W  *  XX(2s6) 

P  *  XX(2s1) 

9  *  XX (2*2) 

R  *  XX(2*3) 

C  DEFINE  THE  TOTAL  LINEAR  ACCELERATIONS  ALONG  THE  MISSILE  XYZ  AXES  AS 
C  UDOTs  VDOTs  WDOT. 

UDOT  *  UMF1  ♦  GAX  -  Q*W  ♦  R*V 

VDOT  *  VMF1  ♦  GAY  -  R*U  ♦  P*W 

WDOT  *  WMFl  ♦  GAZ  -  P*V  ♦  9*U 

XX(1*4)  *  UDOT 
XX(lsS)  *  VDOT 
XX(1s6)  *  WDOT 

C  ZERO  OUT  LINEAR  ACCELERATIONS  ALONG  MISSILE  Y  AND  Z  AXES  UNTIL 

C  ROUND  HAS  CLEARED  LAUNCHER. 

IFCLNEXIT)  60  TO  100 
XX(1*5)  *  0.0 
XX(1s6)  *  0.0 
100  CONTINUE 

C  CALCULATE  AXIAL  AND  TRANSVERSE  MOMENTS  OF  INERTIA  IX*  AND  IY*IZ. 

C  CALCULATE  THE  CENTER  OF  GRAVITY  CG. 

IF(SKIP)  GO  TO  130 
IFC.NOT.  BURNOUT)  GO  TO  125 
IX  *  IXF 
IT  *  IYF 
IZ  *  IY 
IXDOT  «  0.0 
IYDOT  =  0.0 
CG  *  CGF 
MOOT  *  0.0 
XMASS  a  WF/GREFF 
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XMASS  =  XX<2*13) 

SKIP  =  .TRUE. 

GO  TO  130 

125  IX  =  IXO  -  I XDOT *T  IME 
IY  =  IYO  -  I YDOT *TIM£ 

IZ  =  IY 

CG  =  CGO  -  CGDOT*TIME 
130  CONTINUE 

C  CALCULATE  HASS  FLOW  RATE  HOOT  AS  A  FUNCTION  OF  THRUST 

IF(SKIP)  GO  TO  140 
HDOT  =  -<THRSTP)/<ISP*GREFF> 

140  XX(1*13)  =  HOOT 

C  CALCULATE  CHQ* THE  RATE  OAHPING  COEFFICIENT  FOR  RATES  ABOUT  HISSILE 
C  Y  AND  Z  AXES*  AS  A  FUNCTION  OF  HACH  NO. 

CHQ  -  TABLKCCMQT*  HCHQ*  XHACHN*  16) 

C  CALCULATE  THE  OAHPING  MOMENTS (TORQUES)  ABOUT  THE  HISSILE  XYZ  AXES 
C  AS  XHP1 AS*  XHQ1AS*  XHR1AS. 

AUX  =  <DPS*DREF**2)/(2.*VRF) 

XHP1AS  *  P*AUX*CLP 
XH01AS  *  0*AUX*CH0 
XHR1AS  =  R*AUX*CHQ 

C  CALCULATE  CENTER  OF  PRESSURE-CP-AS  FUNCTION  OF  HACH  NO.  AND  ANGLE 
C  OF  ATTACH  ALFTOT. 

CP  *  DREF*TABL2(CPT*  HCP*  CPAF*  XHACHN*  ALFTOT*  8*  2) 

C  CALCULATE  AERODYNAHIC  HOHENT  ARM-AMA 
AHA  -  CP  “  CG 

C  CALCULATE  THE  STABILITY  HARGIN  IN  CALIBERS 
SH  *  AHA/OREF 

C  CALCULATE  THE  AERODYNAHIC  RESTORING  TORQUES  ABOUT  THE  HISSILE  Y  AND 
C  Z  AXES  AS  XHAY*  XHAZ. 

XHAY  *  FW1AP*AHA 
XHAZ  *  -FV1AP*AHA 

C  CALCULATE  TORQUES  ABOUT  HISSILE  XYZ  AXES  DUE  TO  PRODUCT  OF  ANGULAR 
C  VELOCITIES  ABOUT  THE  OTHER  TWO  AXES  TIHES  THE  DIFFERENCE  BETWEEN 
C  INERTIA  OF  THE  OTHER  AXES  AS  TAX#  TAY*  TAZ. 

TAX  *  Q*R*(IZ  -  IY) 

TAY  *  -P*R*< IX-IZ) 

TAZ  *  -P*Q*<IY-IX> 

C  CALCULATE  THE  TORQUES  ABOUT  HISSILE  XYZ  AXES  DUE  TO  TIHE  RATE  OF 
C  CHANGE  OF  INERTIA*  AS  TIX*  TIY*  TIZ.  TORQUE  *  DERIVATIVE  OF 
C  ANGULAR  HOHENTUH  *  (O/DT)  IW  *  I*WDOT  ♦  IDOT*W 
TIX  *  P*IXDOT 
TIY  *  Q*IYDOT 
TIZ  *  R*IYDOT 

C  CALCULATE  THE  TORQUES  ABOUT  THE  HISSILE  XYZ  AXES  DUE  TO  THRUST 
C  COMPONENTS  ALONG  THE  HISSILE  XYZ  AXES*  AS  XMPlGS*  XMQ1GS*  XMRlGS. 
XMP1GS  *  THR$TP*SC ARF 
ARH  *  CG  -  FACTOR 
XHQ1GS  *  ARM*FW1GP 
XMRlGS  »  “ARH*FV1 GP 

C  CALCULATE  THE  TOTAL  TORQUES  ABOUT  THE  HISSILE  XYZ  AXES  AS  TORQX* 

C  TORQY*  TORQZ. 

TORQX  *  XHP1AS  ♦  XMPlGS  ♦  TAX  “TIX 

TORQY  ■  XHQ1AS  ♦  XMQlGS  ♦  XHAY  ♦  TAY  -  TIY 

TORQZ  *  XHR1AS  ♦  XMRlGS  ♦  XHAZ  ♦  TAZ  -  TIZ 

C  ZERO  OUT  THE  TORQUES  ABOUT  HISSILE  Y  AND  Z  ACES  UNTIL  ROUND 

C  HAS  CLEARED  THE  LAUNCHER.  THIS  PUTS  THOSE  ANGULAR  ACCELERATIONS  TO 
C  ZERO. 

IF(LNEXIT)  GO  TO  300 
TORQY  *  0.0 

TORQZ  »  0.0  63 

300  CONTINUE 


C  CALCULATE  ANGULAR  ACCELERATIONS  ABOUT  MISSILE  XYZ  AXES  AS  PDOT, 

C  8D0T ,  AND  ROOT.  THE  ANGULAR  ACCELERATION  ABOUT  ANY  AXIS  IS  THE 
C  TOTAL  T0R8UE  DIVIDED  BY  THE  MOMENENT  OF  INERTIA. 

PDOT  =  TORQX/IX 
8D0T  =  T0R3Y/IY 
ROOT  =  T0R8Z/IZ 

C  XX< 2x1)  IS  THE  ROLL  RATE  AND  WOULD  NORMALLY  BE  CALCULATED  BY 

C  E8UATING  XX<1,1>  TO  PDOT, THE  ANGULAR  ACCELERATION, AND  INTEGRATING 
C  XX<1,1>.  AT  THE  PRESENT  XX<2,1>  IS  CALCULATED  FROM  A  TABLE  OF  ROLL 
C  VERSUS  TIME,  SO  XX <1 ,1 )  WILL  NOT  BE  E8UATED  TO  PDOT. 

XX  <1 ,2)  *  8D0T 
XXC1 ,3)  =  ROOT 

C  CALCULATE  THE  MISSILE  VELOCITIES  IN  THE  GROUND  FRAME~XDOT, YDOT, 

C  ZDOT-AS  XX <1 ,7) ,  XX<1,8>,  XX<1,9> 

DO  305  J=1 ,3 

305  XX(1,J+6)=DCM(1,J)*XX(2,4)+DCM(2,J)*XX(2,5)+DCM(3,J)*XX(2,6) 
PN0M=S8RT<ZMW**2+YMW**2> 

TWOPI  =  2.*PI 
PHII=AMOD(XX<2,20),TWOPI) 

SNPHI=SIN(PHII) 

CSPHI=COS(PHII) 

IF<XMW)310,311,310 

311  IF(PNOM) 312,313,312 

312  ETAPRN=1 .5707963 
GO  TO  314 

313  ET APRN=0 . 0 
GO  TO  314 

310  ETAPRN=ATAN(PNOM/XMW) 

314  ALPHPR=BETPRN*SNPHI+ALPP«N*CSPHI 
BETAPR=BETPRN*CSPHI-ALPPP.N*SNPHI 

C  CALCULATE  THE  EULER  ANGLE  DERIVATIVES  SIDOT,  THDOT ,  AND  FIDOT  AND 

C  EQUATE  THEM  TO  XX(1,10),  XX<1,11),  XX(1,1 2) . 

SIDOT  *  <SFI*Q  +CFI*R)/CTH 
THDOT  =  CFI*Q  -  SFI*R 

FIDOT  =  P  ♦  <SFI*STH*8  ♦  C FI *STH *R> / CTH 

XX  <1 ,1  0)  =  SIDOT 

XX<1,11>  =  THDOT 

XX(1 ,1 2)  =  FIDOT 

JUMP=1 

IFCTIME.GE.TLL. AND. TIME. LT.TEXIT)  50  TO  5030 
IFCTIME.GE.SUSOFF. AND. TIME. LT.TTO)  GO  TO  5030 
IF(K-4>23,8,2060 
3  IFCTIME.LE.O.)  GO  TO  5010 
GO  TO  5020 
5010  NPAG c=Q 
LPAGE^O 

TPR=TIME-TPRINT 
5020  IF ( J UMP. NE.1 )  GO  TO  5030 

IF (TIME- (TPR+TPRINT) >5040,5030, 5030 
5030  ALPTD=RADCON*ALPPRN 
BETTD=RADCON*BETPRN 
ETATD=RADCON*ETAPRN 
ALF= ALPHPR 
3ET0=BETAPR 

IF(XX(2,9).GE.O.)  GO  TO  2002 

TP1*TIME 

XP1=XX(2,7) 

YP1=XX<2,8) 

ZP1=XX<2,9) 

GO  TO  2003  64 

2002  IF(IPRNT.E8.2)G0  TO  2003 
XP2*XX(2,7> 


itlt.-.  ^  J 


YP2=XX(2*6) 

ZP2=XX (2*9) 

TP2=TlME 

IPRNT=2 

DOLL=ZPl  /(ZP2-ZP1) 

XMP=XPl-D0LL*(XP2-XP1) 

YMP=YP1-D0LL*<YP2-YP1) 

TMP=TP1-00LL*<TP2-TP1) 

WRITE(6*2005)  TMP*XMP* YMP 

2005  FORMAT < / *  TIME  OF  IMPACT  = ' *  F3. 4*  '  IMPACT  R AN5E  =  * *  FI  0 . 3* 

1‘  IMPACT  CROSS  RN6=*  *  F9. 3/) 

C  WRITE  THE  VALUES  OF  THE  VARIABLES  AT  GROUND  IMPACT  TO  THE 

C  PLOT  FILE  TAPE2. 

C  THIS  SECTION  OF  CODE  IS  ONLY  REACHED  WHEN  THE  ALTITUDE  XX(2*9) 

C  HAS  GONE  P0STIVE*MEANING  THAT  GROUND  IMPACT  HAS  OCCURED. 

WRITEC2) 

A  TIME*  XX < 2/ 7) *  XX(2*8)*  XX(2*9)*  VRF/ 

8  THRSTP*  XX< 1 *7) *  XX (1  *  8 )*  XX<1*9)*  XMACHN* 

C  XMASS*  XX  <  2*4) *  XX  <2*  5) *  XX(2*6>*  MDOT  * 

D  SM*  XX ( 1  * 4) *  XX<1*5)*  XX(1*6)*  DP* 

6  RHO*  PID*  AID*  RID*  CA* 

F  ALFTOT*  PID1 *  QID1*  RID1*  CN* 

G  ALPTD*  TORCX*  TOR QY *  TORfiZ*  CY* 

H  6ETTD*  WINX*  WINY*  WINZ*  CLP* 

I  X.MQ1AS*  XMAY*  TAY*  TIY*  CMfl* 

J  XMR1AS*  XMAZ*  TAZ*  TIZ*  9  ANK* 

<  SI*  TH,  FI*  UMF1*  VMFl * 

L  WMF1 *  G AX*  GAY*  GAZ*  DIST 

C  SET  THE  LOGICAL  VARIA3LE  IMPACT=. TRUE .  TO  STOp  TRAJECTORY. 

IMPACT  =  .TRUE. 

2003  CONTINUE 

C  CALCULATE  ROLL  RATE  AS  PID  IN  REV/SECS 

PID  =  XX (2*1 ) /TWOPI 
WINX=XX(1*7)+WXYZ(1*1) 

WINY=XXC1*3)-WXYZ<2*1) 

WINZ=WXYZ<3*1)“XX<1*9) 

AL  F0=0 . 

IFCWINX.NE.O.)  ALF0=ATAN2(WINZ*WINX) 

QID=  RADCON*  XX(2*2) 

RID=  R  ADCON* XX ( 2* 3 ) 

PID1  =  RADCON*  XX(1*1) 

DI3T-SQRT((XTAR-XX(2*7)) ** 2+ <Y TAR-X X ( 2 , 8 ) ) ** 2+ ( ZTAR+ XX (2* 9) ) ** 2) 
3ID1=  RADC0N*XX(1*2) 

RID1=  RADCON*  XX<1*3) 

SI  =  XX( 2*10) *RADCON 
TH  =  XX ( 2* 1 1 )*RADCON 
FI  =AMOD(XX(2*12)*TWOPI)*RADCON 
IF (MOD (NpAGE*4) .NE.O)GO  TO  5050 
JR  IT  E ( 6*  5005 ) 

LP4GE=LPAGE+1 
JR  IT  E ( **  5006 )  ID4Y*  LP  AGE 
5050  WP ITE( 6*  5037) 

IFCJJMP. =0.1 )  TPR=TIME 
NPAGE=NPAGE+1 

JR ITt ( 6* 4C00)  TIME*  XX(2*7)*  XX(2*3)*  XX(2*9)*  VRF 
4000  FORMAT (2X** TIME  = • *  El  4. 3* 2X* • R ANGE  = *  *  El  4. 8*2X* • DEFL  ='*E14.3*2> 
' * ' ALT  ='*E14.8*2X*'VELR  ='*E14.8) 

JR ITE ( 6*  401 0  >  THRSTP*  XX(1*7)*  XX(1*3)*  XX(1*9)*  XMACHN 
4010  FCRMAT(2X*’THRUST='*E14.8*2X*,XD0T  = ' *  El  4 . 8  * 2X* • YDOT  =  ’*E14.3*2> 
' * ' ZDOT  ='*E14.8*2X*'MACH  =’*E14.3) 

JRITEC  6*4020)  XMASS*  XX(2*4)*  XX(2*5)*  XX(2*6)*  MDOT 
4020  FORM AT(2X*' MASS  =  • *  El  4. 8* 2X* ’ UUU  =' *E1 4.8*2X* ' VVV  =’*£14.3*2) 
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wTO  C  OvtO  v*  a 


•*'WWW  =' *£14.6*2X* 'MOOT  =  '*E14.3) 

WRITt(6*4030)  $M*  XX(1*4)*  XX(1*5)*  XX(1*6)*  DP 
4030  FORMAT (2X* '  $M  = • * £1 4. 3* 2X* * UDOT  ='*£14.8/2X*'VD0T  =  * 

•*'WDOT  ='*El4.8*2X*'De  =**E14.3> 

WRITE<6*  4040)  RHO,  P ID*  310,  RIO,  CA 
4040  FORMAT (2X* ' RHO  = ' * £1 4 . 3* 2X * • p  =  * /El 4.8*2X** 0  =  ' 

**'R  =**£14.8*2X**CA  =**£14.3) 

WRITE(6*405G)  ALFTOT*  PID1*  QID1*  RID1*  CN 
4050  F0RMAT(2X*'ALFT0T='*£14.8*2X**PID1  = • * E 1 4 . 8* 2X* *  0  ID  1  =  ' 

**'RIDl  ='*E14.8*2X**CN  =**E14.3) 

WRITE ( 6*  4060 )  ALPTD*  TORQX*  TCRQY*  TORCZ*  CY 
4060  F0RMAT(2X*'ALPTD  = '  *  El  4. 3,  2X*  *  TORQX  =  •  *  El  4. 8/ 2X*  '  TORQY  =  • 
'  *  'TORQZ  =**El4.8*2X**CY  ='/E14.8) 

WRITE(6*4C70)  BETTD*  WINX*  WINY*  MINI,  LNEXIT 
4070  FORMAT (2Xr'BETTD  = '  *El  4 . 3* 2X* * WI NX  =  • * El  4 . 8/ 2X* *  WINY  =' 


*S1 4. 8*2> 


=  '*£14. 3*2> 


*  El  4. 3*  2> 


*  El  4.  8*  2> 


=*/E14.8) 

WINY*  MINI,  LNEXIT 


='*E14.8*2X*'WINY  ='*E14.3*2> 


l*El4.8*2X*'LNEXIT='*L3) 


WRITE(6*4030)  XMQ1 AS*  XMAY/  TAY*  TIY*  CMQ 
4080  F0RMAT(2X*'XMQ1A$='*£14.8*2X*'XMAY  =• *E1 4 .8* 2X* ' TAY  =• 

**'TIY  ='*E14.8*2X*'CMQ  =  **E14.3) 

WRITE (6*4090)  XMR1AS*  XMAZ*  TAZ*  TIZ*  BANK 
4090  F0RMAT(2X**XMRlAS='*E14.8*2X**XMAZ  = • *  El  4. 8* 2X* ' TAZ  =' 
1  * ' TIZ  =,*E14.8*2X*'BAN<  =**E14.3) 

WRITE(6*41 00)  SI*  TH*  FI*  UMF1*  VMF 1 
4100  FORMAT (2X**  PSI  =  *  *  El  4.  3*  2X*  '  TH  ET  A  =• *  El  4 . 8* 2 X* ' PHI  = ' 
' * ' UMF1  ='*E14.8*2X*'VMF1  =  **E14.8) 

WR IT  E ( 6*41 1 0 )  WMF1 *  G AX*  GAY*  GAZ*  CGTRVL 
4110  F0RMAT(2X*'WMF1  = ' *  El  4. 3* 2X * ' GAX  = • *  El  4 . 8* 2X* ' GAY  =  ' 

*  * ' GAZ  ='*E14.8*2X*'CGTRVL=**E14. 3) 

C  STOP  PROGRAM  WHEN  IMPACT  IS  TRUE.  IMPACT  IS  TRUE  WHEN 

C  XX (2*9)  GOES  POSTIVE. 

IF(IMPACT) STOP* GROUND  IMPACT’ 

5005  FORMAT <1  HI ) 

5006  FORMAT (5X**  6”D  TRAJECTORY  OF  2.75  MK66I2*  M151  W ARHE  AD ' *  1 
! *12*'/* *12*1  OX* 'PAGE* *12) 

5007  FORMAT (1  HO) 

WRITS<2) 

A  TIME*  XX (2*7) *  XX(2*8)*  XX(2*9>*  VRF* 

B  THRSTP*  XX (1 *7) *  XX(1*8)*  XX(1*9)*  XMACHN* 

C  XMASS*  XX(2*4)*  XX(2*5)*  XX(2*6)*  MOOT* 

D  SM*  XX (1 *4) *  XX(1*5)*  XX(1*6>*  DP* 

E  RHO*  PID*  QID*  RID*  CA* 

F  ALFTOT*  PID1 *  QID1*  RID1*  CN* 

G  ALPTD*  TORQX*  TORQY*  TORQZ*  CY* 

H  BETTD*  WINX*  WINY*  WINZ*  CLP* 

I  XMQ1AS*  XMAY*  TAY*  TIY*  CMQ* 

J  XMR1AS*  XMAZ*  TAZ*  TIZ*  BANK* 

K  SI*  TH*  FI*  UMF1*  VMFl* 

L  WMF1*  GAX*  GAY*  GAZ*  OIST 
C 

C  THESE  EQUATIONS  PERFORM  4TH  ORDER  RUNGS-KUTTA  INTEGRATION 
C 

5040  IF(JUMP-2)23*21*2062 

23  IF (RK  .£Q.  2.0)  SO  TO  300 
IF (K  .EQ.  4)  KN  =  0 
KN=KN+1 
DO  22  J  =1 *N 

CK(KN*J)=XX(1*J)*DELTT 
22  CONTINUE 

GO  TO  <407*409*411 *420) *K 
420  K=1 

DO  404  1=1*2  66 

DO  404  J=1 *N 


*  El  4 . 3  *•  2> 


*  El  4. 3*  2) 


*  £1  4. 8*  2> 


=  '  *  El  4.  3*  2> 


ALTITUDE 


0X*I2*'/' 


404  Y(I,J)=XX(I,J) 

00  405  J  =NN,  N 

405  XX(2,J)=CK(1,J)/2.+Y(2,J) 

TIME=TIME+D£L02 

50  TO  413 
407  K=  2 

00  408  J=NN, N 

403  XX(2,J)=CK(2,J)/2.+Y(2,J) 

GO  TO  41 3 

409  <=  3 

DO  410  J  =NN, N 

410  XX(2,J)=CK(3,J)*Y(2,J) 

TIME=TIME+0EL02 

50  TO  413 

411  <=4 

00  412  J  =NN, N 

412  XX  (2, J)=<CK<1, J)+2.*CK<2,J)+2.*CK<3,J)*CK(4,J) )/6.+Y(2,J) 

413  CONTINUE 
GO  TO  598 

C  THESE  EQUATIONS  PERFORM  2N0  ORDER  RUNGE-KUTTA  INTEGRATION 

800  CONTINUE 

GO  TO(9CO,  850)  <2 
850  K2  =  1 

00  860  1=1/2 
00  860  J  =  NN,  N 
360  YCI/J)  =  XX(I,J) 

C  Y(1,J)  =  DERIVATIVE  AT  THE  BEGINNING  CF  THE  TINE  STEP  AND 

C  Y ( 2, J ) =THE  STATE. 

00  870  J  =  NN,  N 

870  XX<2,J)  =  Y ( 2, J  )  +  0.6666657*0GLTT*Y(1,J) 

TIME  =  TIME  ♦  DELTT*0. 6666667 
GO  TO  1000 
900  K2  =  2 

00  910  J  =  NN,  N 

910  XX(2,J)  -  Y(2,J)  +  0.25*OELTT*U(1,J)  ♦  3.*XX(1,J)) 

TIME  =  TIME  +  DELTT*0. 3333333 
GO  TO  1000 
593  CONTINUE 

C  CALCULATE  LOGICAL  VARIABLE  BURNOUT 

BURNOUT  =  TIME  .GE.  SUSCFF 
1000  IF (HEIGHT)  600,  1010,  1010 
1010  IF (RK  .EO.  4.0)  GO  TO  1020 
IFCTIME  -  TFINAL)85,85,700 
1020  I F (TIME  -  TFINAL)7,7,700 

700  HR  IT  E (6, 701 ) 

701  FORMAT  ( *  TIME  EXCEEDS  MAX, RUN  NEXT  CASE') 

600  J UMP-2 

C  TIME  OF  FLIGHT  HAS  EXCEEDED  ALLOTED  RUNTIME  TFINAL.  SET  NOTIME 

C  EQUAL  TO  TRUE  TO  STOP  PROGRAM. 

NOTIME  =  .TRUE. 

GO  TO  8 

21  STOP*  FLT .TIME  GREATER  THAN  RUN  TIME* 

ENO 

SPROG  TABLK 
SXREF 

FUNCTION  TABLK(X,Y,XNUMB, LENGTH) 

DIMENSION  X< 30) , Y< 30) 

00  1  1=2, LENGTH 

IF(XNUMB-Y(I))3,3,1  67 

1  CONTINUE 
I=LENGTH 

3  TA8LK*X(I-1)+<XNUM3-Y<I-1))*<X<I)-X<I-1))/(Y(I)-Y(I-1)) 


SPROG 

SXREF 


SPROG 

SXREF 


RETURN 

END 

MULTY 

SUBROUTINE  MULT Y ( A/B/N/ C ) 

DIMENSION  A(3/3)/B(3/3)/C<3/3) 

DO  1  L  =  1/N 
DO  1  1=1/3 
C(I/L)=0. 

DO  2  L  =  1/N 
DO  2  1=1/3 
DO  2  J=1/3 

C(I/L)=C<I/L)*A(I/J)*B(J/L> 

RETURN 

END 

INVERT 

SUBROUTINE  INVERTCA/C) 

DIMENSION  A( 3/3) /C (3/ 3) 

CC1/1)=A(2/2)*A(3/3)-A(2/3)*A(3/2) 

C(1/2)=A(1/3)*A<3/2)-A<1/2>*A(3/3) 

C(1/3)=A(1/2)*A(2/3)-A<1/3)*A<2/2) 

C(2/1)=A(2/3)*A<3/1>-A(2/1)*A<3/3) 

C(2/2)=A(1/1)*A<3/3)-A<1/3)*A(3/1) 

C(2/3)=A(1/3)*A(2/1)“A<1/1)*A<2/3) 

CC3/1)=A<2/1)*A(3/2)-AC2/2)*AC3/1) 

C(3/2)=A(1/2)*A(3/1)-AC1/1>*A(3/2) 

C(3/3)=A(1/1)*A(2/2)“AC1/2)*A(2/1) 

DETER=A(1/1)*C(1/1)+A<1/2)*CC2/1)+AC1/3)*C(3,1) 

DO  5  1=1/3 

DO  5  J=1/3 

C(I/J)=C(I/J)/DETER 

RETURN 

END 

TA3L2 

FUNCTION  TABL2(TAB/TMAC/TALF/X/Y/N01/N02) 

DIMENSION  TA9<17/5)/TMAC(30>/TALF(30) 

DO  1  I*2/N01 

IFCX. LT.TMACCD)  GO  TO  2 
CONTINUE 

I  =  N01 
CONTINUE 
DO  3  J=2/N02 

IFCY. LT.TALF(J))  60  TO  4 
CONTINUE 

J  =  N02 
CONTINUE 

PALF=(Y-TALF(J-1 )>/<TALF(J)-TALF(J-1 )) 
PMAC=(X-TMAC(I-1))/<TMAC<I>-TMAC<I-1)) 
P1=TAB(I-1/J-1)+PALF*(TAB(I-1/J)-TABCI-1/J-1 )) 

P2=TAB(I/J-1 )+PALF*<TA3(I/J)-TAB(I/J-1 )) 

TABL2=P1+(P2-P1)*PMAC 

RETURN 

END 

HELVEL 

SUBROUTINE  HELVELCXC ALL/ YC ALL/ZC  ALL/VX/ VY/VZ/TIME/ IREAD) 

COMMON  / KEEP/  RAD/ XNR8/0RL/NUMREC/RR0T/DT/ PERIOD 
COMMON/FUSE/  NPTS/  VXINF<160)/  VZINF(160) 

I/X8ARC160)/  Y8AR<160)/  Z3AR<160)/  SIGXC160)/  SIGZC160) 

5/  DDJ (160/4) /  ET AK(1 60/4) /  XIK(l60/4) 

68 


.*•  -v  >v  -N  -N  -N  -■»  •  *• 


SPROG 
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3/ 

XLX<160)/ 

XMX<1 60)/ 

XNX  <1 60) 

4/ 

XLE<160)/ 

XME  <1 60) / 

XNE  <1 60) 

5/ 

XLZ<1 60) / 

XMZC160)/ 

XNZ  <1 60) 

COMMON  /VORT/  X < 6/90/ 2 ) / Y< 6/ 90/ 2 ) / Z < 6/ 90/2 ) / A <6/ 90/ 2 ) / 

1  GAMA<0/90/2)/SEG<6/90/2)/3AMB1 <100)/VXF/VYF/VZF/NPS<6)/ 

2  NA/NB/NW/NA3/XMCL/XMSL 

DIMENSION  DUM <90/2) /DUMP <160)/ DUMMY <1 00) /DUMY< 6) 

IF (IREAD  .NE.  0)  GO  TO  50 

RAD=. 01745329252 

REWIND  10 

READ (1 0/11 5)  NREC 

NUMREC=NR£C+1 

REWIND  10 

DO  10  1  =  1 / NUMR  EC 

RE AD  <1 0/ 1 1 5)  IDUM 

RE AD <10/155)  <<X(I/J/K)/J=1/90)/K=1/2)/<(Y(I/J/K)/J=1/90)/K=1/2) 
*/<<Z(I/J/IO/J=1/90)/K=1/2>/<<GAMA(I/J/K)/J=1/90)/K=1/2) 
*/<<A<I/J/K)/J=1/9G)/K=1/2)/<<SEG<I/J/K)/J=1/ 90) /K=1 / 2) 

READ <10/ 135)  PSI/PSIF/DP3I/NPS<I) 

)  CONTINUE 

READ <10/1 45)  D/D/RR0T/D/D/NR9/NPTS 
RE AD  <1 0/1 55)  OUMM/DUMM/D/D/D/D 
READ<1 0/165)  D/NA/NB/NAB/NW/D 

READ <10/1 55)  D/D/D/XMCL/XMSL/XLAM/D/D/D/D/D/D/D/D/D/DUMMY/GAMB1/ 


115 

135 

145 

155 

165 


SPROG 

SXREF 


READ <10/1 45)  D/D/RR0T/D/D/NR9/NPTS 
RE AD  <1 0/1 55)  OUMM/DUMM/D/O/D/D 
READ<1 0/165)  D/NA/NB/NAB/NW/D 

READ <10/1 55)  D/D/D/XMCL/XMSL/XLAM/D/D/D/D/D/D/D/D/D/DUMMY/GAMB1/ 
*XBAR/Y9AR/ZBAR/SIGX/SIGZ/DDJ/XIK/ETAK/XLX/XMX/XNX/XLE/XM£/XNE/ 
*XLZ/XMZ/XNZ/ DUMP Y/DUMMY/DUMY/VXINF/VZ INF 
FORMAT  <11 0) 

i  F0RMAT<3<1E12.5)/I10) 

F0RMAT<5<1PE12.5)/2I8) 

F0RMAT<6<1PE12.5)) 

F0RMAT<1PE12.5/4I10/1PE12. 5) 

XNRB=NRB 

0RL=XNRB*6.*RAD*RR0T*XLAM 
REWIND  10 

PE RI0D=60. /XNR3/ FLOAT <NB) 

DT=PSRI0D/ FLOAT <NUMREC) 

RTIME=AMOD<TIME/ PERIOD) 

IF  <RTINE  .LT.  0.)  RTIME=RTIME+PERIOD 

DEX=RTIME/DT 

FRAC*AM0D<DEX/1 .) 

FR AC  C=1 • ~FR AC 

M=DEX+1 

XX=XCALL/RR0T 

YY=YCALL/RROT 

ZZ-ZCALL/RROT 

CALL  FUSLGE<XX/YY/ZZ/VXF/VYF/VZF) 

CALL  VLCTY<XX/YY/ZZ/V1/V2/V3/M) 

CALL  VLCTY<XX/YY/ZZ/V11/V22/V33/M*1) 

VX=<V1*FRACC+V11 *FRAC)*ORL 

VY=<V2*FRACC+V22*FRAC)*0RL 

VZ=<V3*FRACC+V33*FRAC)*0RL 

RETURN 

END 

5  VLCTY 

SUBROUTINE  VLCTY < XIPT/ YIPT/ Z IPT/VX/VY/VZ /MM) 

COMMON  /VORT/  X < 6/90/ 2 ) / Y < 6/ 90/ 2 ) / Z < 6/ 90/2 ) / A< 6/ 90/ 2) / 

1  6AMA<6/90/2)/SEG<6/90/2)/GAMBl<100)/VXF/VYF/VZF/NPS<6)/ 

2  NA/NB/NW/NAB/XMCL/XMSL 
M=M0D<MM-1/6)+1 

VX=0.  69 

VY  =  0. 


-AOi.  v  vV  v.  A  v.v.'.  v  v  v  v.v  .  -.v. 


W’J 


^V.T.TT.jr.T-. 


vz=o. 

XXX=XIPT 
YYY=YIPT 
ZZZ=ZIPT 
00  SO  J=1/N3 

SIG2=SflRT<<XXX-X01/1/J>>**2+(YYY-Y<M/1/J>>**2'*'(ZZZ-ZC!1/1/J>>**2> 

00  40  K=1 /NW 

SIG1=SIG2 

SIG2*SQRT<(XXX-X<M/K+1/J)>** 2+ (YYY-Y (M/K+l / J ) ) **2 
*  +  <ZZZ-Z<M/K-M/J>>**2> 

SEGSa=SEG<M/</J>**2 
HM1  =  SIG1**2+SIG2**2 
I F (HM1 . GT . S  E  G  S  3) GO  TO  30 

HM2  =  .25*<(SIG1*SIG2)**2-SEGSa)*(SEoSa-(SIG1-SIG2)**2)/SEGSa 
I F  C  Hr12  .GT.  A(M/K/  J  )  **2)  GO  TO  30 
GGG=GAMA(M/K/J)/S£G(M/</J> 

GO  TO  31 

GGG=GA?1A(M/K/J)*<SIGl+SIG2)/{SIGl*SIG2*((SIGl  +  SIG2)**2“SEGSQ)) 
DIFX=XXX-X<M/K/J> 

DIFY=YYY-Y(M/K/J) 

DIFZ=ZZZ-Z<M/K/J) 

XNU1*DIFY*<Z<M/K/J)-Z<N/K  +  1/J>>-DIFZ*(YCM/K/J>-YCM/KM/J>> 

XNU2  =  DIFZ*<X<M/K/J)-XOl/K*1/.)>>-DIFX*<Z<M/IC/J>-2<M/K  +  1/J>> 

XNU3  =  0IFX*<Y<M/K/J)-Y<M/K  +  1/J>)-DIFY*CX  ("!/</  J > “X  CM/ K*1 / J > > 

VX  =  VX  +  XNU1 *GGG 

VY  *  V Y  «-XNU2*GGG 

VZ  =  VZ  +XNU3*GGG 

CONTINUE 
CONTINUE 

SIG1  *  SQRTCXXX**2+YYY**2+ZZZ**2> 

00  60  L*1 / N3 

LPS=M0D<NPS<M)+CNA*<L-1)>/NB*NAB/NA> 

IF (LPS.EQ.O)  LPS  *  NA 

XNU1*ZZZ*Y<M/1/L)"YYY*Z(M/1/L) 

XNU2=XXX*Z<M/1/L)-ZZZ*X<M/1,L> 

XNU3=YYY*X<M/1/L)-XXX*YCM/1/L> 

SIG2  =  SQRT<<XXX-XCM/1/L)>**2+<YYY-YCM/1/L)>**2*CZZZ-ZCM/1/l>>W> 


OEN  =  SIG1*SIG2*<<SIG1+SIG2)**2-1.0) 
IF(OEN.EQ.O.O)  OEN  >  .001 
GGG  =  GAMB1(LPS)*CSIG1+SIG2)/DEN 
VX=VX+XNUl*GGG 
VY=VY+XNU2*GGG 
VZ=VZ+XNU3*GGG 
60  CONTINUE 

VX=VX  +  XMCL  +VXF 
VY*VY  ♦VYF 
VZ*VZ  -XMSL  +VZF 
RETURN 
END 

SPROG  FUSLGE 
SXREF 

SUBROUTINE  FUSLGE<XD/YD/ZD/UFD/VFD/WFD> 
COMMON/FUSE/  NPTS/  VXINFC160)/  VZINF<160> 

1 / X9 A  R ( 1 60) /  YBAR(1 60)/  ZBARC160)/  SIGXC160)/ 

2/  OOJ (160/4)/  ETAKC1 60/4) /  XIK(l60/4) 

S/  XLXC160)/  XMXC160)/  XNX(160> 

6/  XLE( 1 60) /  XME (1 60)/  XNE(160> 

7/  XLZC160)/  XMZC160)/  XNZ(160) 

DIMENSION  R J ( 4) /  EJ<4)/  HJ(4>/  EMJC4) 

SUMU  *  0.0 

SUMV  *  0.0  70 

SUMW  =  0.0 


SIGZC160) 


GO  TO  13 


IF (NPTS. £9.0) 

DO  12  J  =1 x  NPTS 
NFL6  =  1 
XB  =  XD-XBAR(J) 

T9  =  YD-Y3AR  < J ) 

ZB  =  ZD-ZBAR(J) 
06=<XIK<J*4)-XIK<J*2))**2+<ETAK<J*4>-ETAK<J*2>>**2 
D5=<XIK<J*3)-XIIUJ*1>>**2+<£TAK<J*3>-£TAK<J*1>>**2 
D7  *  AMAXl <D5*D6> 


XI 

=  XLX(J) 

* 

XB 

♦ 

XMX(J) 

*  Y  8 

♦ 

XNX(J) 

*  ZB 

ETA 

=  XLE(J) 

* 

XB 

+ 

XME(J) 

*  Y3 

♦ 

XNE(J) 

*  Z  3 

ZETA 

*  XLZ ( J  ) 

* 

XB 

♦ 

XMZ< J) 

*  Y8 

♦ 

XNZ(J) 

*  Z  3 

RO  =  XI**2+ETA**2+ZETA**2 
TJ  =  R0/D7 

IF(PO.NE.O.O)  GO  TO  4 

VXI  =  0.0 

VETA  =  0.0 

VZETA  =  6.2831853072 

GO  TO  10 

IF  (  TJ.LT.6.)  GO  TO  5 

SJ=.5*<XIK<J*3>-XIK<J*1>>*<£TAK<J*2>-ETAK<J*4>>/<R0*S0RT<R0>> 

VXI  =  SJ  *  XI 

VETA  *  SJ  *  ETA 

VZETA  *  SJ  *  ZETA 

GO  TO  10 

VXI  *  0.0 

VETA  =  0.0 

VZETA  *  0.0 

DO  9  K*1*2 

DO  9  1*1*4 

II  «  I  ♦  1 

IF(I.Ea.4)  II  *  1 

TRMX  =  XIK ( J  *  1 1 )  ”  XIK<J*I) 

TRHE  *  ETAK<J*I1)-ETAK<J*I> 

IFOC.EQ.2)  GO  TO  7 

TRM1  *  XI  -  X  IK ( J  *  I ) 

TRM2  *  ETA  -  ETAK<J*I> 

H J ( I )  *  TRM1  *  TRM2 

EJ (I )  *  TRM1 **2  ♦  ZET A**2 

R J (1 )  *  <TRM2**2  ♦  EJ  < I) ) ** • 5 

XF(TRMX.EQ.Q.G)  TRMX  =  1.0E-6 

EMJ  <  I) *  TRUE  /  TRMX 

GO  TO  9 

TRMR  *  (RJ(I)+RJ(I1)-DDJ(J*I))/(RJ(I)+RJ(I1)+DDJ(J*I)) 

IF(TRMR.LT .0. )  WRITE <6*9 9)  J* I*XD*YD* ZD*TRMR*RJ < I) *R J < II  > *DD J  <J*I> 
FORMAT  <  4X*52HTRMR  IS  NE  G— J  *  I  *X  D*  YD*Z  D*  TRMR*  R  J  <  I)  *R  J  <  II  >  *DD  J  <  J  *  I. 
1  / <21 5*7E1 5.5)) 


TRMR 

3 

ALOG (ABS (TRMR)  ) 

TRME 

s 

TRUE 

/ 

DDJ  <J*I) 

TRMX 

=-TRMX 

/ 

DDJ <J*I) 

VXI 

3 

VXI 

TRUE  *  TRMR 

VETA 

3 

VETA 

TRMX  *  TRMR 

I F (ZET A. EO . 0 .0)  GO  TO  9 

VZETA  *  VZETA  ♦  ATAN< (EMJ (I) *EJ  < I )  -HJ(I)  ) / (ZET  A*R  J  < I)  )> 
*  ATAN((£NJ<I)*EJ(I1)-HJ<I1))/<ZETA*RJ(I1))) 

CONTINUE 

VVX  *  XLX(J>*VXI*XLE(J>*VETA*XLZ<J>*VZETA 
VVY  =  XMX<J)*VXI*XME<J>*VETA*XMZ(J)*VZ£TA 
VVZ  *  XNX(J)*VXI+XNE(jj*VETA+XNZ(J)*VZETA 
IFCNFLG.EQ.2)  GO  TO  11 

VVVX  *  VVX 
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v.v. 


vvvz  =  vvz 

YB  =  -Y0-Y3AR<J) 

NFLG  =  2 
GO  TO  3 

11  TRM  =  SIGXCJ)*VXINF<J>+SIGZ<J>*VZINFCJ> 
SUMU  =  SUMU+TRM*<VVVX+VVX) 

SUMV  =  SUMV+TRM*<VVVY-VVY> 

SUMU  *  SUMW+TRM*(VVVZ+VVZ) 

12  CONTINUE 

13  UFO  =  SUMU 
VFD  =  SUMV 
WFO  *  SUMU 
RETURN 

END 

SBEND 
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APPENDIX  E 
RPLOT  Listing 


•••••••  THIS  program  writes  data  from  unformatted  plot  file 

~DIMENSION  A (60) 

REAL  MO 


DATA  YES,NO/3HYESt3H  NO/ 
DPUTIun  lK 


CALL  CONNEC (5) 

_ CALL  CONNECi&X _ 

SO  FORMAT (5X**IS  OATA  TAPE  FORMATTED  ?  ENTER  YES  OR  NO*) 

UOTTr / A .CA 


51  FORMAT <A3) 


IF (ANS  .FQ.  NO)GO  TO  70 

_ READ  (ISLfciRfXS _ 

46  FORMAT (l0X**NUMB£R  OF  RECORDS 


52  FORMAT(l0X.*VA®lAdLE  NUMRfR*,5X  **MINIMUM*,13X**MaXIMUM*) 

—  WRITE  (20,523 _ 

WRITE (6*100) 

_ READ  L5  *2001 MV _ 

47  FORMAT (10X*I6»l4X*ElA»8*feX»E14.8) 


READ (15) xmjn,XMAX 

_ W^IT£120,421IjXMINiAMAX _ 

60  CONTINUE 

_ GO  10. 75  _ _ _ _ _ 

70  CONTINUE 

— lDQ-EQRMJLri5jL**fUI££  2...UIG1I  NUMtiER=NU.flF_I)ARlABLJES  QN-EACM.  .RE 
WRITE (6*100) 

_.20JQ_  FQFMATIIZI _ 

RE AO (5*200)  NY 

_ 75  DQNILNUE _ 

300  FORMAT <5X**ENTER  TST0P* F6.3*TO  STOP  READING  PLOT  TAPE  AT*) 


WRITE (6*300) 

_ WRITE  (6*301) _ 

400  FORMAT (F6» 3) 

_ REAI1I5*.400)_TSJDR 

425  FORMAT ( lHl ) 


t.q 

R 

5  READ ( 1 5) (A(N),K=1»NV) 

_  IFiEOF  il5J  )60*-L0 _ 

10  WRITE  (20» 5000  (A(i\)*K=l»Nv) 


♦3X6HZ7Z  =,E15.8*3X6HVEL  R^El^.S/TH  MASS  =*E15.8*3X6HXD0T 
*E15.  8  »  3X6HYDDT  =»£L5. 9.*5A6HZP0I.  =  *£15^8i6A6HMA£H_=uE15«aZ- 
*7H  THRST=,E15.R*3X6HOPP  =*E15.8,3X6HQ0U  =*E15.8,3X6HPRR 

*  El  5-  6  ♦  3  X6HQT  N?.R=  »E  L5^87.Z  H  .PI  ST  *El5-.J9j3A6HgQDI  J?*£15*EL* _ 

*3X6HQf)UT  =,E15.9,3X6HRDDT  =»E15.8,3X6HALFI  =*E1&.8/7H  CG 
*E15,8.3X6HUUU  =*E15.8,3X6HVVV  =,E15.8*3X6HWWW  =*E15.A* 


1.^  Os”.’ 


«  '  4  *. » 


APPENDIX  F 


TEKPLOT  Listing 


PREVIOUS  PAGE 
IS  BLANK 


IIW 


*•  »'j- ••  v v v 


PROGRAM  TEKPLOT(TAPE2«TAPf:S=6E>»TAPE6=65* TAPE  10 ♦OUTPUT) 


DATA  YES,NO/3HYESt3H  NO/ 

Al  I 


CALL  RINITT 


CALL  ERASE 


CALL  RELL 

UT-3 


1  REWIND  NT 


CALL  HOME 


►****THE  FOLLOWING  STATEMENT  SETS  THE  RUN  NUMBER  TU  ONE* 


The  NUMBER  of  VARIABLES  ON  THF  DATA  RUN.  IF  MORE  THAN  TWENTY*  CHANGE 


130  FORMAT (12) 


I  E 


140  FORMAT (2E14.8) 


STANDARD  FORMAT. 


WRITE(6*131> 


WRITE (6*133) 


WRITE (6*134) 


RR 


,%■ 

v 

WRITE (6*137) 

|v 

WRITE(6*138> 

M 

H im  liiKuii 

r#' 

k 

WRITE (6*139) 

% 

139  FORMAT (5x**IF 

\%vvVr*.s’CsCNV.v.-C..-..  s  •  •  •  • 


145  FORMAT (A3) 

READ  (5*14S)  IANS1  -  -  - —  - 

IF(IANS1.EQ.YES)GU  TO  5 

- CALL-  FORMAT-Qqy  *NT4-N£IS3 - - 

DISPLAY  A  MESSAGE  ASKING  T He.  OPERATOR  TO  ENTER  A  2  DIGIT  NO.  =  TO  THE  NO 

-OF— P-LQT S  TU.  BE  MADE _ - _ 

5  CALL  HOMF 

~CALL  RELl^ 

_ CALL  -AMMOCg _ 

WRITE  <6*131 ) 

- WflXXEt^»4S0) - - — 

mo  FORMAT(l0X.*ENTER  a  2  DIGIT  NO.=TO  NO.  DF  PLOTS  TO  BE  MADE*) 

-RFAO  THE.  _QP£RATURS  ANSWER - 

READ (5» 160) NPLTS 

— lAO-  FORMAT  t-12) _ 

CALL  ERASE 

_ CALL  HOME _ 


TO  BE  PLOTTED 


CALL  BELL 
WpITE 

by  ENTERING  PAIRS  OF  2  DIGIT  NUMBERS  CORRESPONDING  TO  THE  INDEPENDENT 

-AND-  OEPENDEMT^-VARIABLE^-PALrS - I - 

CALL  ANMODE 


WRITS (6.13) ) 


WRITE (6.161) 

.FQRMAT-CUX.ffJ 
WRITE (6» 162) 

4.62! 

WRITE (6*163) 

163  FORMAT  tlflX.*  OF  THF  PATH  nFSTGNATFB  THF  iNn.  VARIARLF*) 
WRITE (6*164) 

WRITE (6*165) 


WRITE (6*166) 

-L66.5jQRMATI10D(.*THF  NUMBER  Of  .ELUIS— IH-BE.  .M&DE.  QR  YOU  Will  RE  STUCK  ») 
WRITE (6*167) 

JLfi7_EQRMA.TjUL0A«*IN— A  RFAQ-LUOP-*! _ 

READ  THE  OPERATORS  ANSWER  FROM  THE  SCREEN 

_ 00  20  I  =  1.NPLTS _ 

?0  RE AO (5* 170) INDV ( I ) * IDEPV ( I ) 

170  FORMAT  (I?«1X«T?) _ 

PROMPT  OPERATOR  TO  ENTER  TITLE  FOR  PLOTS 

_ CALI _ HOME _ 

CALL  ERASE 


_ CALL  RELI _ 

CALL  ANMODE 

171  FORMAT  tlx. »FNTFR  TITLE  FOP  PI  OTS-HP  TO  80  CHARACTERS* ) 
WRITE (6* 171 ) 

_ _  R£Afl(5tl46) TITLF _ 

146  FORMAT(BAIO) 

BEGIN__PLQT-T-LNG _ 

1*0  CONTINUE 

_ 00-500  .  Ial  .NPI  T.S _ 


*>  ttt:  t-  r 


•*  <-  • 


.T^v  .  Tr  ,TT,  ?■,  V,  .’.  ; 


■  V  W  V  ~J-  ".> ' J?  w. 


IF (IANS »EQ* YES) GO 
CAU  RELL _ 


TO  S 


CALL 

CALL 


HOME 
ANMQOE 


WRITE (6» 131 ) 
WRITF (ftigSQ) 


250  FORMAT (20X,*PLOTTING  COMPLETE-*) 
_ CAU.  FTNITT (0.400) _ 


STOP 

END 


SUBROUTINE  XYGRID (XMIN,XMAX» YMIN,YMAX> 
■(imE.M5LUN-XXlLa.YX3.) _ 


X(l)=2. 

jmx=2« 


X(2)=0. 

Y_L21=XU_ 


X  (3)  =0  • 
Y(3)=ft. 


*  v  v 


*3 


>'A 


v  .  ' 


vV*j 

4  ^  V 


— * 


«• 


CALL  ANMQOE 


/  / 


>v*  s  ■*»_!;  TV'S.*'  ■>  \i  *m  v  ■*  >j  » j  .mini  i  i_«p  y  7 ?  r».Tyy 

to 

I 


REWIND  NT 

_ .READ  t  NI )  NPI5 _ 

DO  190  KK=1*NV 

TOO  -REAP  (Ml)  XMI-N-(KK)  *XMAX  WU - - 

SELECT  THE  SUBSCRIPTS  OF  THE  INDEPENDENT  AND  DEPENDENT  VARIABLE 

_  JUsINOM(I) _ 

K2=IDEPV ( I ) 

-Call  subroutine  XYfiiaajHis  draws,  the  or  id  and_lebels  ijl.wjtm- ihe-mih» 

VALUES  OF  the  INDEPENDENT  * DEPENDENT  VARIABLE 

_ CALL  F1RASE _ 

CALL  HOME 

THE  FIRST  VARIABLE  ON  A  RECORD* A(l) *IS  USUALLY  TIME.  IF  YOU  WANT 
THE  LOWEST  VALUE  ON  THE  HORIZONTAL  AXIS  TO  BE  ZERO* REGARDLESS  OF  THE 

-SMALLEST  VALUE  OF  T lME*-ACTj  V  ATE- THE— FOLLOWING  -STATEMENTS - 

IF(KI.NE.I)  GO  TO  200 

- —  - - 

XMA=XMAX (Kl)-XMlN(Kl) 

- CALL— XVfiRlO  (XMI-«XMA*XMIN  (K21  »XMAX  (KZJ  ) - 

*  GO  TO  201 

-200  CON-TI-NUE -  -  -  -  —  —  - - 

CALL  XYGRID(XMIN(M)  *XMAX(KI)  *XMIN(K2)  ,XMAX(K2) > 

-201  -CONTJNUE - 

MOVE  THE  CURSOR  TO  THE  STARTING  POINT  OF  THE  GRAPH*A(K1> * A(K2) 

DRAW  THF  -SREC1F-1E 0-JGftAP.B - 1 - 

DO  220  lA=l*NPTS 

- READ (NT) ( A  ( -J ) ,Js.l»NV) - 

IF (EOF (NT)) 221*222 

222  CONTINUE - 

IF  THE  FIRST  VARIARLE*A(1) *USUALLY  TIME*00ES  NOT  START  AT  ZERO 
-YQILCAN— INSFPT-A  STATFmFNT  THaL  SURRA CTS  OFF  THE.  MINIMUM  -VALUE- OF- XI ME 
AND  THE  ZERO  TIME  ON  THE  PLOTS  WILL  REPRESENT  THE  STARTING  TEST  TIME 

_ SUCH  A  STATEMENT— WQUI  D  BF — IF  (Kl  *EQ.  1 )  A  (K1 ) -A  (Kl ) -XMLMU1- _ 

IF(IA.EQ.1)CALL  MOVEA ( A  (K] ) ,A(K2)) 

_ TF(IA.FQ.l)GO  TQ-22Q _ i _ 

CALL  DRAWA (A (K1)*A(K2>) 

-220  CONTINUE - 2 - 

221  CONTINUE 

LABEL  THF -X  AND- Y  . AXIS - 

CALL  LABLR(K1*K2*TITLE) 

- - CALL-MOCOPY - - 

CALL  BELL 

SO a  CONTINUE _ 

DISPLAY  A  MESSAGE  ON  THE  TEKTRONIX  SCREEN  ASKING  THE  OPERATOR  IF  HE 

-WANTJSl.  TO- CONTINUE  PLOTTING  FRqm.  THE-CURRENT  RUN _ 

CALL  HOME 

CALL  ERASE _ _ 

CALL  ANMODE 
WRITE (6i?3Q) 

230  FORMAT (10X**00  YOU  WISH  TO  CONTINUE  PLOTTING  FROM  THE  CURRENT  RUN* 

_ U - 

WRITE (6*240) 

-2A0  -F_flRMAT  (lQX*» ENTER  YES-OR  -NOJL) — - - 

READ  THE  OPERATORS  ANSWER  FROM  SCREEN 

_  ftF  AQ^SalASHAMS--..  - -  _  _  — 
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-  a.-  , 


'  *  %  > 


.'*■  »'* 


‘a'.Va’ 


m 


m 


REWIND  2 


9(1) =10000000. 

— s..cm=-iooooooo. _ 

6  READ (2) (A(J) *J=1*NV) 

- IF-1EQF-12) )  2Q»  1 0 _ 

10  NPTS=NPTS*1 
FIND  THE  MlN.*B(I) .AND  THE  -MAX*  «C  ( I ) *0F  EACH  VARIABLE 
DO  15  1=1 *NV 

- atn=AHlNl  (AlDtBlI)) - 

15  C(I)=AMAXl(A(I).C(I)) 

_ on  to  ft _ 

20  REWIND  2 

- NTs  10 _ : - 

WRITE (NT )NPTS 

- 00  .30  1  =  1  *NW - 

30  WRITE (NT) B ( I ) *C ( I ) 

- QO  60  I?1«NPTS _ 

READ (2) (A ( J) *J=1*NV) 

Afl-WRIIE  (NT UAI J1  »  J=-l  >  N  V ) _ 

FNDFILE  NT 

_ REMIND  NT _ 

RETURN 

_ END _ 

BLOCK  DATA 

_ COMMON  /LARI /t.AR(60*3) _ 

DATA  , (LAB(I*J) *J=1*3) *1=1*19) 

— A/10HTIME  OE-EL»HmifihI  IN  SF.1QHCS. - 

B10H  RANGE  *10H  *10H 

_ C10HDFF1,ECIIUN»10H-PUSTT VF  R.  1 0 HIGH. I - 

010HALTITUDE  .ldH  *10H 

_ E10HRELATIVE-.V»1QHELUCITY _ tlOH _ 

F10HMA9S  IN  SL.10HUGS  *10H 

_ G10HRANGE— VELQ«lflHCLTY _ «1QH _ 

HIOHCROSS  RANG* 10HE  VEL0ClTY*10H  POS-RIGHT 
JUOHVERTICAL  V*10HELOCITY-PO*10HSTIVE  DOWN 
J10HMACH  NUMBE* 10HR  -*10H 

— K10KTHRUST - »iOH - - 

L10HROLL  RATE  *10HIN  REV/SEC*10H 

— M10HP-ITCrt-RATE*UH-IN  DEGS/S*10HEC - 

N10HYAW  RATE  I.10HN  OEGS/SEC*10H 

_ OlOHOYNAKIC  PR.  1  flHFSSHRE _ »10H - - — 

P10HRANGE  TO  G.IOHO  TO  TARGE*10HT 

— 01 QHROLI.  ACCEL  »1&H»  IN  RAD/StlOHFC/SFC - 

R10HPITCH  ACCE* 10HL.  IN  RAD/*10HSEC/SEC 

— S10HYAM  ACCF(.*»  lt)K--lN  RaQ/SF* 1 OHC/SFC - 1— 

DATA  ( (LAB (I*J)*J=1*3) *  1=20*38) 

_ A/10HALFI _ _ U-OB - %. 

B10HCENTER  OF  *10HGRAVITY-FT*10H.  W/R  TAIL* 

— CIOHMISSILE  AX.10H1S  UELUClT.l QHl-U - *— 

D10HMISSILE  AX.10HIS  VELUClT*lOHY-V  * 

riAUMTCCTir  A V .  1 AUI C  UPl  UCrT.lAuV.U  . 


F10H0ENSITY  *10H 


*  10H 


DATA  ( (LAB (I * J) * J=1 »3) » 1=19*57) 


A10HCL 


•  10HH 
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